Bioinformatics Analysis of Hub Genes and Potential Therapeutic Agents Associated with Gastric Cancer.

PURPOSE: The current treatment methods available for advanced gastric cancer are not very promising. Hence, it is important to explore novel biomarkers and potential therapeutic agents to treat gastric cancer (GC). This study aimed to identify hub genes associated with GC prognosis and explore potential drugs for its treatment.
MATERIALS AND METHODS: Three gene expression data of GC and normal tissues were downloaded from the Gene Expression Omnibus (GEO) and processed to identify the differentially expressed genes (DEGs). We conducted a comprehensive analysis of DEGs, including functional enrichment analysis, construction of protein-protein interaction (PPI) network, identification of hub genes, survival analysis and expression verification of hub genes. Finally, we constructed the network of miRNA-mRNA, and predicted the drugs that might be effective for GC treatment.
RESULTS: A total of 340 DEGs, including 94 up-regulated and 246 down-regulated genes, were identified. Among the up-regulated DEGs, the enrichment terms were primarily related to tumorigenesis and tumor progression, extracellular matrix organization, and collagen catabolic process. Additionally, 10 hub genes (FN1, COL3A1, COL1A2, BGN, THBS2, COL5A2, THBS1, COL5A1, SPARC, and COL4A1) were identified, out of which 7 genes were significantly associated with poor overall survival (OS) in GC. The expression levels of these 7 hub genes were verified using real-time PCR, immunohistochemistry, and the GEPIA2 (Gene Expression Profiling Interactive Analysis) server. A regulatory network of miRNA-mRNA was also constructed, and the top 4 interactive miRNAs (hsa-miR-29b-3p, hsa-miR-140-3p, hsa-miR-29a-3p, and hsa-miR-29c-3p) that targeted the most hub genes were identified. Finally, fourteen small molecules were predicted to be effective in treating GC.
CONCLUSION: The identification of the hub genes, miRNA-mRNA network, and potential candidate drugs associated with GC provides new insights into the molecular mechanisms and treatment of GC.

Introduction

Gastric cancer (GC) is one of the leading causes of death worldwide,1–3 with the highest incidence in East Asia, followed by eastern Europe.4 Early gastric cancer lacks specific symptoms, which leads to the fact that many gastric cancer patients are already in advanced stage when diagnosed.5 Heterogeneity of GC is also one of the factors leading to the low survival rate of patients with GC.6 Molecular diagnosis and targeted therapy are effective means to improve the early diagnosis rate of gastric cancer,7 and molecular diagnosis must be based on the search for specific tumor biomarkers.8

With the completion of the Human Genome Project, cancer has increasingly been studied at the genetic level. DNA microarrays play an important role in identifying of genes leading to cancer and have the characteristics of high throughput, accuracy, sensitivity, and low cost. The technology is already being widely used in disease diagnosis and drug screening.9,10

In this study, three gene expression profiles of GC and normal tissues were downloaded from the Gene Expression Omnibus (GEO) and processed to identify the differentially expressed genes (DEGs). Then we conducted a comprehensive analysis of DEGs, including functional enrichment analysis, construction of protein–protein interaction (PPI) network, identification of hub genes, survival analysis and expression verification of hub genes. Finally, we constructed the network of miRNA-hub genes, and predicted the drugs that might be effective for GC treatment through the online database.

Materials and Methods

Acquisition of Data

Three microarray datasets (GSE54129, GSE79973, and GSE13911) containing GC and normal tissues were downloaded from GEO online database (, accessed time: April 24, 2021) using GEOquery and preprocessed using affy R packages. All microarray datasets based on GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) contained more than 10 samples.

Screening of DEGs

The probe set IDs were converted into gene symbols using the annotation package hgu133plus2.db (3.2.3). If multiple probe sets corresponded to the same gene, the mean value of those probe sets expression was used. R software (version 4.0.3) was used to process all data using the Limma (3.46.0) package and Biobase (2.50.0).11,12 Fold change was considered significant for P<0.05 and threshold values >1.0 or <−1.0.

GO and KEGG Enrichment Analysis of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID) (, version 6.8, accessed time: April 28, 2021) provides a comprehensive set of functional annotation tools for investigators to understand the biological significance of a list of genes. The results of functional enrichment analysis were visualized by by R software package ggplot2.13

PPI Network Construction and Module Analysis

PPI network analysis was performed using the online tool Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, , accessed time: May 11, 2021), and visualized using Cytoscape software (version 3.8.2). Then several major modules of the network were analyzed using MCODE plug-in of Cytoscape software, and hub genes were identified using Cytoscape plugin cytoHubba.

Survival Analysis and Online Validation of Hub Gene Expression

The Kaplan–Meier plots of the top 10 hub genes were drawn using the Kaplan-Meier plotter (, accessed time: May 14, 2021), and the survival prognosis forest map was visualized using the randomForest package.

The expression levels of hub genes between GC and normal tissues were compared using Gene Expression Profiling Interactive Analysis (GEPIA2; , accessed time: May 15, 2021). Furthermore, the genetic mutation information for hub genes was explored using the online tool of cBio Cancer Genomics Portal (, accessed time: May 16, 2021).

Real-Time PCR

We collected GC and para-cancer tissues from thirteen patients who had been first diagnosed with GC at the Endoscopy Center of the First Affiliated Hospital of Guangxi Medical University from July to November 2021. None of these patients had received any tumor-related treatment or had any other malignancy. Total RNA was isolated using TRIzol (R0016, Beyotime). Complementary DNA (cDNA) was synthesized using PrimeScript™ RT reagent Kit (RR047A, Takara). The qPCR reaction was set up using SYBR green PCR master mix containing 0.5 µL of cDNA. The quantitative real-time PCR analysis was performed using the real-time PCR system ABI 7500 Fast (Applied Biosystems). Data were normalized to GAPDH, and fold change in mRNA expression relative to controls was calculated using the 2−ΔΔCt method.

Immunohistochemistry (IHC)

We retrospectively studied the postoperative tumor tissues and normal tissues from 60 patients diagnosed with GC in The First People’s Hospital of Suqian from 2019 to 2020. The patients did not have any other known tumors, and no radiotherapy or chemotherapy had been performed before surgery. After removing paraffin, hydration, and sealing, the specimens were mixed with the antibodies anti-FN1 (#26836, Cell Signaling Technology), anti-COL1A2 (ab270994, Abcam), anti-BGN (ab209234, Abcam), anti-THBS2 (ab112543, Abcam), anti-COL5A1 (D264482, Sangong Biotech), anti-SPARC (ab225716, Abcam), and anti-COL4A1 (D260202, Sangong Biotech), and incubated overnight at 4 °C (with dilution ratios of 1:500, 1:2000, 1:2000, 1:1000, 1:200, 1:100, and 1:100, respectively). DAB (3, 3ʹ-diaminobenzidine) was used for signal detection. Additionally, we performed statistical analysis of images from immunohistology studies using Image-Pro Plus 6.0 software. Finally, we use GraphPad Prism 8 to visualize the statistical results.

Diagnostic Value of Hub Genes Based on TCGA Data

We downloaded the RNA sequencing dataset of stomach adenocarcinoma (STAD) from The Cancer Genome Atlas (TCGA, , accessed time: May 17, 2021). We compared the expression of hub genes in tumor tissues and non-tumor tissues of GC patients to evaluate their diagnostic value in GC, and the ROC curve was visualized using the online tool Helix (, accessed time: May 18, 2021).

Construction of miRNAs-mRNA Network

The Encyclopedia of RNA Interactomes (ENCORI, ; version 3.0, accessed time: May 21, 2021) is a public platform that identifies more than 2.5 million miRNA–mRNA interactions from multi-dimensional sequencing data. The target miRNAs of hub genes were screened based on the following criteria: SLIP-DATA ≥ 5, Degradome-Data ≥ 0, pan-Cancer ≥ 6, and programNum ≥ 3. Additionally, the miRNAs-hub genes network was visualized using Cytoscape software (version 3.8.2).

Drug-Hub Gene Interaction

The DGIdb (, accessed time: May 25, 2021) contains over 40,000 genes and 10,000 drugs involved in over 100,000 drug-gene interactions and belonging to one of the 42 potentially druggable gene categories. The drugs with potential therapeutic effect on the hub genes of GC were screened by DGIdb. Additionally, we used the online tool STITCH (, accessed time: May 25, 2021) to construct an interaction network between hub genes and potential drugs.

Statistical Analyses

We conducted the paired samples t-test to analyze mRNA expression and immunohistological positive rate of hub genes in tumor tissues and paracancerous non-tumor gastric tissues with IBM SPSS Statistics v.24 software (IBM, Chicago, IL, USA). A value of P <0.05 was regarded as statistically significant.

Results

The details of the three datasets downloaded from GEO are presented in Table 1. As shown in Figure 1, a total of 340 DEGs were identified in the three microarray datasets, of which 246 were down-regulated (Figure 1A) and 94 up-regulated (Figure 1B). The gene names of these 340 DEGs are shown in Table 2.

Table 1

Detailed Information on the GEO Microarray Profiles of GC Patients

GEO Profile	Type	Source	Case	Control	Platform	Annotation Platform
GSE54129	mRNA	Gastric cancer	111	21	GPL 570	Affymetrix Human Genome U133 Plus 2.0 Array
GSE79973	mRNA	Gastric cancer	10	10	GPL 570	Affymetrix Human Genome U133 Plus 2.0 Array
GSE13911	mRNA	Gastric cancer	38	31	GPL 570	Affymetrix Human Genome U133 Plus 2.0 Array

Abbreviation: GEO, Gene Expression Omnibus.

Figure 1

The 340 DEGs identified from the three microarray datasets: (A) 246 downregulated genes; (B) 94 up-regulated genes.

Table 2

Details of the 340 DEGs

DEGs	Gene Names
Up-regulated (94)	RCN3, GREM1, LAMA5, COL3A1, WASF1, OLFML2B, APOC1, TGM2, AMIGO2, COL4A2, DUSP14, SRPX2, HOXC9, PLPPR4, COL5A1, FXYD5, SULF2, LEF1, LGALS1, FAM83D, BAG2, NNMT, THY1, HOXB6, CDH11, KLHL23, CEBPB, CPXM1, GAL, CST4, DNM1, COL4A1, COL1A2, FJX1, CDH13, THBS1, ITGA5, PI15, TIMP1, SPHK1, CLEC5A, TNFAIP6, SNX10, SPOCK1, SPARC, HOXC10, CST2, PLA2G7, COL5A2, HOXB7, PRRX1, IL13RA2, PLAU, NT5DC2, FCGR1B, COL6A3, IGF2BP3, CRISPLD1, ASPN, AJUBA, NID2, RARRES1, TREM2, PMEPA1, COL12A1, SFRP4, NOX4, BGN, MMP11, TMEM158, CTHRC1, CLEC11A, FN1, HOXC6, MFAP2, SULF1, THBS2, FKBP10, CLDN1, CCN4, COL11A1, CXCL8, LY6E, FAP, COL8A1, SPP1, ADAMTS2, SERPINH1, CEMIP, COL10A1, CDH3, CST1, FNDC1, INHBA
Downregulated (246)	ATP4B, CBLIF, ATP4A, AQP4, GKN1, LIPF, ESRRG, SOSTDC1, LTF, KCNE2, TRIM50, GKN2, KCNJ16, CHGA, MFSD4A, CWH43, CCKBR, PGC, TMED6, MAP7D2, GC, CXCL17, MUC5AC, VSIG2, MUCL3, FBP2, FAM3B, CLIC6, KRT20, VSIG1, ADH1C, SCNN1B, AKR1B10, RFX6, IGH, RDH12, CPA2, SLC26A9, DNER, FCGBP, ALDH3A1, ABCA8, MAL, SULT2A1, KCNJ15, SSTR1, HPGD, GUCA2B, HHIP, TNFRSF17, CNTN3, NKX2-3, LINC02381, ADAM28, LIFR, ELAPOR1, ADTRP, CAPN9, B3GNT6, SLC28A2, CA2, AADAC, ADH7, AKR1C2, PIK3C2G, PSAPL1, GATA6-AS1, GCNT2, HAPLN4, FBXL13, PKIB, ALDOB, ADGRG2, CYP2C18, AKR1C1, SMIM5, SULT1C2, HEPACAM2, CAPN13, UBE2QL1, SMIM24, LYPD6B, SCIN, RNASE1, SPINK7, AXDND1, VSTM2A, MGAM, LDHD, MST1L, RPS6KA6, ENPP6, SPINK2, UGT2B15, ZBTB7C, PAIP2B, UNC5B-AS1, PBLD, CYP3A5, ATP2A3, SPTSSB, CYP2C9, BCAS1, CA9, B4GALNT3, CAPN8, AMPD1, AKR1C3, INSC, PRDM16, JCHAIN, TRPA1, GCKR, METTL7A, ALDH1A1, CLTRN, ITPKA, MRAP2, LINC01133, VILL, ASPA, NRG4, PDGFD, ACKR4, ADRB2, APOBEC1, GPAT3, RAB27A, PROM2, TMEM220, SMIM6, AKR7A3, TENT5C, TREH, LRRC31, TMEM238L, SLC14A1, FMO5, SLC9A1, CYP3A7-CYP3A51P, TMPRSS2, TMEM38A, SELENBP1, ARHGEF37, SLC41A2, ADAMTS15, PLAAT2, FER1L4, NQO1, CSTA, RAB37, SULT1B1, MYRF, RAB27B, SYTL2, PPP1R36, UPK1B, DGKD, NTN4, SIDT2, RASSF6, CYP4F12, ATP13A4, RNASE4, ADH1A, FA2H, OASL, DDX60, KBTBD12, SLC7A8, CCL28, SH3BGRL2, CYSTM1, STYK1, NOSTRIN, TMEM266, ACE2, ANKRD36BP2, PWWP3B, LRRC66, GRAMD1C, MAOA, TMEM171, EPB41L4B, PPFIBP2, NKX2-2, XK, ANG, SALL1, FAM174B, ETFDH, STS, MIOS-DT, PCAT18, STX19, KLF4, S100P, IQGAP2, HGD, EPN3, PER3, PLLP, OSBPL7, CLDN23, CIDEC, DHRS7, NR3C2, PELI2, CYP2C19, GOLGA2P10, ABCC5, PLCXD1, C1orf116, DISP1, TMEM116, GALE, LNX1, UBL3, PTPRH, SMPD3, GGT6, C4orf19, PAQR8, KLHDC7A, ARHGEF28, COBLL1, ENTPD5, HAPLN1, GPRC5C, RNF128, TEX9, DAPK2, SLC22A23, CDHR2, GPD1L, KLB, MAGI3, KCNK10, ABCA5, TMEM92, LOC101928881, CYP3A43, FBP1, PXMP2, PTGR1, EIF2AK3

To better understand the biological functions of overlapping DEGs, their enriched biological processes (BP), cellular components (CC), and molecular functions (MF) were studied. For the BP category, the up-regulated DEGs were involved in the extracellular matrix organization, collagen catabolic process, collagen fibril organization, cell adhesion, and endodermal cell differentiation (Figure 2A), whereas downregulated DEGs were significantly involved in steroid metabolic process, xenobiotic metabolic process, oxidation-reduction process, digestion, and daunorubicin metabolic process (Figure 2E). For the CC category, the up-regulated DEGs correlated with the extracellular region, extracellular matrix, extracellular space, proteinaceous extracellular matrix, and endoplasmic reticulum lumen (Figure 2B), whereas the downregulated DEGs were associated with extracellular exosome, organelle membrane, integral component of membrane, apical plasma membrane, and basolateral plasma membrane (Figure 2F). For the MF category, the up-regulated DEGs were enriched in extracellular matrix structural constituent, extracellular matrix binding, integrin binding, platelet-derived growth factor binding, and calcium ion binding (Figure 2C), whereas the downregulated DEGs were related to oxidoreductase activity, retinol dehydrogenase activity, monooxygenase activity, oxygen binding, and aromatase activity (Figure 2G). For the KEGG pathway, the up-regulated DEGs were enriched in ECM-receptor interaction, focal adhesion, protein digestion and absorption, and the amoebiasis PI3K-Akt signaling pathway (Figure 2D), and the down-regulated DEGs were enriched in chemical carcinogenesis, metabolism of xenobiotics by cytochrome P450, retinol metabolism, drug metabolism-cytochrome P450, and gastric acid secretion (Figure 2H). The detail results of GO and KEGG pathway analysis of up-regulated and down-regulated DEGs are shown in Tables 3 and 4, respectively.

Figure 2

The top 10 significant terms in GO analysis (Biological Process, Cellular Components, and Molecular Function) and KEGG enrichment analysis of DEGs. (A–D) Down-regulated DEGs; (E–H) Up-regulated DEGs.

Table 3

GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) Enrichment Analysis Results of Up-Regulated DEGs

Category	Term	Count	P value	Genes
GO_BP	Extracellular matrix organization	20	4.17E-19	LAMA5, OLFML2B, SPARC, COL11A1, BGN, FN1, NID2, THBS1, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, MFAP2, SPP1, COL6A3, COL8A1, COL10A1, ITGA5
GO_BP	Collagen catabolic process	13	5.87E-16	COL11A1, COL12A1, MMP11, COL3A1, ADAMTS2, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, COL6A3, COL8A1, COL10A1
GO_BP	Collagen fibril organization	10	4.06E-13	GREM1, MMP11, COL3A1, ADAMTS2, COL1A2, COL5A1, COL11A1, COL12A1, COL5A2, SERPINH1
GO_BP	Cell adhesion	17	2.12E-09	TNFAIP6, COL12A1, FN1, THY1, THBS2, NID2, THBS1, CDH3, COL5A1, FAP, CDH11, SPP1, SPOCK1, COL6A3, COL8A1, CDH13, ITGA5
GO_BP	Endodermal cell differentiation	7	5.06E-09	COL4A2, COL11A1, COL12A1, FN1, COL8A1, INHBA, ITGA5
GO_CC	Extracellular region	37	1.27E-15	LAMA5, CEMIP, SPARC, OLFML2B, CXCL8, COL11A1, COL12A1, TREM2, THBS2, PLA2G7, NID2, THBS1, ADAMTS2, PLAU, SPP1, COL10A1, TIMP1, IL13RA2, CLEC11A, FN1, BGN, INHBA, SFRP4, MMP11, COL3A1, GAL, COL1A2, COL4A2, COL5A1, FNDC1, CRISPLD1, COL4A1, COL5A2, MFAP2, APOC1, COL6A3, COL8A1
GO_CC	Extracellular matrix	19	4.93E-15	LAMA5, COL12A1, BGN, FN1, THBS2, NID2, THBS1, ASPN, MMP11, COL3A1, LGALS1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, COL6A3, COL8A1, TGM2
GO_CC	Extracellular space	32	1.11E-13	LAMA5, SPARC, CXCL8, TNFAIP6, COL12A1, THBS1, CST4, CST2, CST1, SRPX2, LGALS1, PLAU, SERPINH1, SPP1, SPOCK1, CPXM1, TIMP1, IL13RA2, CTHRC1, CLEC11A, FN1, SULF1, SULF2, GREM1, SFRP4, COL3A1, GAL, COL1A2, FAP, COL6A3, CDH13, FJX1
GO_CC	Proteinaceous extracellular matrix	17	2.48E-13	OLFML2B, SPARC, COL11A1, BGN, FN1, ASPN, MMP11, ADAMTS2, LGALS1, COL1A2, COL5A1, COL5A2, SPOCK1, COL6A3, COL10A1, TIMP1, CTHRC1
GO_CC	Endoplasmic reticulum lumen	15	5.78E-13	FKBP10, COL11A1, COL12A1, THBS1, RCN3, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, SERPINH1, COL6A3, COL8A1, COL10A1
GO_MF	Extracellular matrix structural constituent	8	4.72E-08	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL5A2, BGN
GO_MF	Extracellular matrix binding	6	1.84E-07	SPARC, OLFML2B, COL11A1, SPP1, BGN, THBS1
GO_MF	Integrin binding	8	1.06E-06	LAMA5, COL3A1, COL5A1, FAP, FN1, THY1, ITGA5, THBS1
GO_MF	Platelet-derived growth factor binding	4	2.05E-05	COL3A1, COL1A2, COL5A1, COL4A1
GO_MF	Calcium ion binding	14	6.53E-05	FKBP10, SPARC, SULF1, THBS2, NID2, THBS1, ASPN, SULF2, RCN3, MMP11, CDH3, CDH11, SPOCK1, CDH13
KEGG	ECM-receptor interaction	14	3.44E-16	LAMA5, COL11A1, FN1, THBS2, THBS1, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, SPP1, COL6A3, ITGA5
KEGG	Focal adhesion	14	2.99E-11	LAMA5, COL11A1, FN1, THBS2, THBS1, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, SPP1, COL6A3, ITGA5
KEGG	Protein digestion and absorption	10	7.28E-10	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL12A1, COL5A2, COL10A1, COL6A3
KEGG	Amoebiasis	10	3.91E-09	LAMA5, COL3A1, CXCL8, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL5A2, FN1
KEGG	PI3K-Akt signaling pathway	14	1.76E-08	LAMA5, COL11A1, FN1, THBS2, THBS1, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, SPP1, COL6A3, ITGA5

Table 4

GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) Enrichment Analysis Results of Down-Regulated DEGs

Category	Term	Count	P value	Genes
GO_BP	Steroid metabolic process	10	1.40E-09	SULT1B1, CYP2C9, CYP3A7-CYP3A51P, AKR1B10, UGT2B15, AKR1C3, AKR1C2, CYP2C19, CYP3A5, SULT2A1
GO_BP	Xenobiotic metabolic process	12	1.72E-09	NQO1, SULT1B1, CYP2C9, ALDH3A1, CYP3A7-CYP3A51P, AKR7A3, AADAC, AKR1C1, UGT2B15, CYP2C19, CYP3A5, CYP2C18
GO_BP	Oxidation-reduction process	26	2.99E-08	HPGD, MAOA, HHIP, ADH7, CYP2C19, CYP2C18, LIPF, AKR7A3, RDH12, LDHD, NQO1, HGD, AKR1C1, CYP3A43, AKR1C3, AKR1C2, PTGR1, CYP4F12, FMO5, CYP2C9, ALDH3A1, FA2H, DHRS7, AKR1B10, ALDH1A1, GPD1L
GO_BP	Digestion	10	4.80E-08	CAPN8, CAPN9, AKR1B10, CCKBR, AKR1C1, GKN1, AKR1C2, SSTR1, PGC, SULT2A1
GO_BP	Daunorubicin metabolic process	4	8.53E-05	AKR1B10, AKR1C1, AKR1C3, AKR1C2
GO_CC	Extracellular exosome	67	3.73E-08	CLIC6, FCGBP, TREH, JCHAIN, SULT1C2, UPK1B, UBL3, AKR7A3, CA2, FAM3B, SH3BGRL2, TMEM38A, ENPP6, BCAS1, HGD, PLLP, ENTPD5, LIFR, RNASE4, PTGR1, RNASE1, PGC, MUC5AC, EPN3, ACE2, C1ORF116, SLC7A8, SCNN1B, CDHR2, ALDH1A1, ANG, ALDOB, FBP1, ASPA, GPRC5C, FBP2, LTF, SMIM24, CSTA, HPGD, TMPRSS2, IQGAP2, SLC9A1, PBLD, SELENBP1, ADGRG2, SCIN, COBLL1, PDGFD, METTL7A, GC, PROM2, NQO1, MGAM, GGT6, AKR1C1, AKR1C3, GUCA2B, RAB27A, RAB27B, GALE, AKR1B10, S100P, CYSTM1, GPD1L, SLC26A9, CCL28
GO_CC	Organelle membrane	11	9.94E-08	FA2H, CYP2C9, CYP3A7-CYP3A51P, AADAC, CYP3A43, UGT2B15, CYP2C19, CYP4F12, CYP3A5, CYP2C18, FMO5
GO_CC	Integral component of membrane	93	1.38E-05	KLB, HEPACAM2, ATP2A3, KLHDC7A, AQP4, CWH43, TMEM92, TMEM220, SMPD3, TMEM266, UPK1B, STS, VSIG1, FAM174B, DNER, CLDN23, CA9, TMEM38A, TMED6, PAQR8, CYP3A7-CYP3A51P, KCNK10, PLLP, DAPK2, ABCC5, CYP3A43, UGT2B15, LIFR, TMEM171, DISP1, SLC22A23, ACE2, ALDH3A1, SLC7A8, RNF128, CCKBR, B3GNT6, NRG4, SCNN1B, MAL, ADTRP, VSTM2A, GPRC5C, SIDT2, GRAMD1C, SMIM24, SMIM5, KCNE2, MYRF, MAOA, SLC41A2, STX19, PELI2, ETFDH, TMPRSS2, ADRB2, SMIM6, CYP3A5, SLC9A1, MFSD4A, ADAM28, ADGRG2, GCNT2, TNFRSF17, SPTSSB, GPAT3, METTL7A, B4GALNT3, PROM2, SLC14A1, NQO1, MGAM, GGT6, FER1L4, TRPA1, ABCA5, STYK1, KCNJ15, KCNJ16, ABCA8, LRRC66, FMO5, CYP4F12, FA2H, ATP4A, XK, AADAC, TMEM116, PXMP2, ATP13A4, CYSTM1, SLC26A9, MRAP2
GO_CC	Apical plasma membrane	13	0.000236875	MGAM, RAB27A, RAB27B, ADRB2, CYP4F12, SLC9A1, UPK1B, ADGRG2, SCNN1B, CDHR2, MAL, SLC26A9, PROM2
GO_CC	Basolateral plasma membrane	10	0.000359953	SLC14A1, HPGD, SLC7A8, CA2, KCNJ16, CA9, AQP4, DISP1, SLC9A1, PROM2
GO_MF	Oxidoreductase activity	15	8.51E-08	HPGD, ADH1C, MAOA, AKR1C1, AKR1C3, ETFDH, AKR1C2, ADH7, PTGR1, CYP2C19, CYP3A5, CYP2C9, ALDH3A1, DHRS7, ALDH1A1
GO_MF	Retinol dehydrogenase activity	5	4.65E-05	RDH12, ADH1C, ADH1A, AKR1C3, ADH7
GO_MF	Monooxygenase activity	7	5.40E-05	CYP2C9, CYP3A7-CYP3A51P, CYP3A43, CYP2C19, CYP4F12, CYP3A5, CYP2C18
GO_MF	Oxygen binding	6	0.000202147	CYP2C9, CYP3A7-CYP3A51P, CYP2C19, CYP4F12, CYP3A5, CYP2C18
GO_MF	Aromatase activity	5	0.000246051	CYP3A7-CYP3A51P, CYP3A43, CYP4F12, CYP3A5, CYP2C18
KEGG	Chemical carcinogenesis	12	9.25E-09	CYP2C9, ALDH3A1, CYP3A7-CYP3A51P, ADH1C, ADH1A, CYP3A43, UGT2B15, ADH7, CYP2C19, CYP3A5, CYP2C18, SULT2A1
KEGG	Metabolism of xenobiotics by cytochrome P450	11	5.41E-08	CYP2C9, ALDH3A1, AKR7A3, ADH1C, ADH1A, AKR1C1, UGT2B15, AKR1C2, ADH7, CYP3A5, SULT2A1
KEGG	Retinol metabolism	10	1.83E-07	CYP2C9, CYP3A7-CYP3A51P, RDH12, ADH1C, ADH1A, UGT2B15, ALDH1A1, ADH7, CYP3A5, CYP2C18
KEGG	Drug metabolism-cytochrome P450	10	3.14E-07	CYP2C9, ALDH3A1, ADH1C, MAOA, ADH1A, UGT2B15, ADH7, CYP2C19, CYP3A5, FMO5
KEGG	Gastric acid secretion	9	6.32E-06	ATP4B, ATP4A, KCNE2, CCKBR, KCNK10, CA2, KCNJ15, KCNJ16, SLC9A1

PPI Network and Module Analysis

The PPI network contained 321 nodes and 664 edges, and its PPI enrichment P-value was lower than 1.0e-16 (Figure 3). In addition, an interaction network of 340 DEGs and their neighboring genes was also constructed (Figure 4A). According to the ranking of MCODE scores, we selected the top three modules for visualization (Figure 4B, module 1, MCODE score=16.111; Figure 4C, module 2, MCODE score=4.167; and Figure 4D, module 3, MCODE score=4). As shown in Figure 4B, module 1 consisted of 19 nodes and 145 edges. KEGG enrichment analysis showed that genes in module 1 were mainly enriched in ECM-receptor interaction, focal adhesion, protein digestion and absorption, PI3K-Akt signaling pathway, and amoebiasis pathway (Table 5). Module 2 consisted of 13 nodes and 25 edges (Figure 4C), with genes enriched in retinol metabolism, chemical carcinogenesis, metabolism of xenobiotics by cytochrome P450, steroid hormone biosynthesis, and drug metabolism-cytochrome P450 pathway (Table 5). Module 3 contained 4 nodes and 6 edges (Figure 4D), with no gene enriched in any pathway (Table 5). The PPI enrichment P-value of each module was lower than 0.05.

Figure 3

A PPI network of the 340 DEGs in GC. The red represents up-regulated genes, and the blue represents down-regulated genes.

Figure 4

(A) Interaction network of the 340 DEGs and their neighboring genes. Red represents selected genes, and greed represents interacting genes. (B–D) The top 3 significant clusters of modules.

Table 5

GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) Enrichment Analysis of Each Module

Category	Term	Count	P value	Genes
Module 1
GO_BP	Extracellular matrix organization	16	4.67E-27	SPARC, COL11A1, BGN, FN1, NID2, THBS1, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, COL6A3, COL8A1, COL10A1, ITGA5
GO_BP	Collagen catabolic process	11	1.32E-20	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL12A1, COL5A2, COL10A1, COL6A3, COL8A1
GO_BP	Collagen fibril organization	7	1.91E-12	COL3A1, COL1A2, COL5A1, COL11A1, COL12A1, COL5A2, SERPINH1
GO_BP	Endodermal cell differentiation	6	6.13E-11	COL4A2, COL11A1, COL12A1, FN1, COL8A1, ITGA5
GO_BP	Cell adhesion	9	1.01E-08	COL5A1, COL12A1, FN1, COL6A3, COL8A1, ITGA5, THBS2, THBS1, NID2
GO_CC	Endoplasmic reticulum lumen	13	2.31E-20	COL11A1, COL12A1, THBS1, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, SERPINH1, COL6A3, COL8A1, COL10A1
GO_CC	Extracellular matrix	14	3.33E-20	COL12A1, BGN, FN1, THBS2, NID2, THBS1, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, COL6A3, COL8A1
GO_CC	Extracellular region	17	1.64E-15	SPARC, COL11A1, COL12A1, BGN, FN1, THBS2, NID2, THBS1, COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL5A2, COL6A3, COL8A1, COL10A1
GO_CC	Collagen trimer	9	1.30E-14	COL3A1, COL1A2, COL5A1, COL11A1, COL12A1, COL5A2, SERPINH1, COL10A1, COL6A3
GO_CC	Proteinaceous extracellular matrix	9	7.59E-11	SPARC, COL1A2, COL5A1, COL11A1, COL5A2, BGN, FN1, COL10A1, COL6A3
GO_MF	Extracellular matrix structural constituent	8	3.45E-13	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL5A2, BGN
GO_MF	Platelet-derived growth factor binding	4	1.67E-07	COL3A1, COL1A2, COL5A1, COL4A1
GO_MF	Extracellular matrix binding	4	2.61E-06	SPARC, COL11A1, BGN, THBS1
GO_MF	Integrin binding	5	4.04E-06	COL3A1, COL5A1, FN1, ITGA5, THBS1
GO_MF	Collagen binding	4	3.35E-05	SPARC, SERPINH1, FN1, NID2
KEGG	ECM-receptor interaction	12	5.27E-20	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL5A2, FN1, COL6A3, ITGA5, THBS2, THBS1
KEGG	Focal adhesion	12	9.83E-16	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL5A2, FN1, COL6A3, ITGA5, THBS2, THBS1
KEGG	Protein digestion and absorption	10	4.12E-15	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL12A1, COL5A2, COL10A1, COL6A3
KEGG	PI3K-Akt signaling pathway	12	3.09E-13	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL5A2, FN1, COL6A3, ITGA5, THBS2, THBS1
KEGG	Amoebiasis	8	2.69E-10	COL3A1, COL1A2, COL4A2, COL5A1, COL4A1, COL11A1, COL5A2, FN1
Module 2
GO_BP	Xenobiotic metabolic process	7	3.74E-12	CYP2C9, ALDH3A1, AKR1C1, UGT2B15, CYP3A5, CYP2C18, CYP3A7
GO_BP	Steroid metabolic process	4	2.54E-06	CYP2C9, UGT2B15, CYP3A5, CYP3A7
GO_BP	Oxidation-reduction process	5	0.000414041	CYP2C9, ALDH3A1, AKR1C1, ALDH1A1, CYP2C18
GO_BP	Drug catabolic process	2	0.003924595	CYP2C9, CYP3A5
GO_BP	Lipid hydroxylation	2	0.003924595	CYP3A5, CYP3A7
GO_CC	Organelle membrane	5	1.56E-07	CYP2C9, UGT2B15, CYP3A5, CYP2C18, CYP3A7
GO_CC	Endoplasmic reticulum membrane	5	0.001255792	CYP2C9, UGT2B15, CYP3A5, CYP2C18, CYP3A7
GO_CC	Intracellular membrane-bounded organelle	4	0.003920426	CYP2C9, ALDH3A1, UGT2B15, CYP3A5
GO_CC	Extracellular space	4	0.042447579	ALDH3A1, CXCL8, PLAU, ANG
GO_MF	Oxygen binding	4	3.29E-06	CYP2C9, CYP3A5, CYP2C18, CYP3A7
GO_MF	Oxidoreductase activity	5	5.91E-06	CYP2C9, ALDH3A1, AKR1C1, ALDH1A1, CYP3A5
GO_MF	Monooxygenase activity	4	6.23E-06	CYP2C9, CYP3A5, CYP2C18, CYP3A7
GO_MF	Heme binding	4	8.23E-05	CYP2C9, CYP3A5, CYP2C18, CYP3A7
GO_MF	Iron ion binding	4	0.000114202	CYP2C9, CYP3A5, CYP2C18, CYP3A7
KEGG	Retinol metabolism	6	1.45E-08	CYP2C9, UGT2B15, ALDH1A1, CYP3A5, CYP2C18, CYP3A7
KEGG	Chemical carcinogenesis	6	4.52E-08	CYP2C9, ALDH3A1, UGT2B15, CYP3A5, CYP2C18, CYP3A7
KEGG	Metabolism of xenobiotics by cytochrome P450	5	2.47E-06	CYP2C9, ALDH3A1, AKR1C1, UGT2B15, CYP3A5
KEGG	Steroid hormone biosynthesis	4	6.55E-05	AKR1C1, UGT2B15, CYP3A5, CYP3A7
KEGG	Drug metabolism-cytochrome P450	4	0.000105515	CYP2C9, ALDH3A1, UGT2B15, CYP3A5
KEGG	Metabolic pathways	7	0.003353172	CYP2C9, ALDH3A1, UGT2B15, ALDH1A1, CYP3A5, CYP2C18, CYP3A7
Module 3
GO_BP	Anterior/posterior pattern specification	4	1.04E-07	HOXC9, HOXB7, HOXB6, HOXC6
GO_BP	Embryonic skeletal system morphogenesis	3	1.57E-05	HOXC9, HOXB7, HOXB6
GO_BP	Transcription, DNA-templated	4	0.001575955	HOXC9, HOXB7, HOXB6, HOXC6
GO_BP	Regulation of transcription, DNA-templated	3	0.022617461	HOXC9, HOXB7, HOXB6
GO_CC	Nucleus	4	0.026223769	HOXC9, HOXB7, HOXB6, HOXC6
GO_MF	Sequence-specific DNA binding	4	2.87E-05	HOXC9, HOXB7, HOXB6, HOXC6
GO_MF	Transcription factor activity, sequence-specific DNA binding	3	0.009344922	HOXB7, HOXB6, HOXC6

Ten up-regulated genes (FN1, COL3A1, COL1A2, BGN, THBS2, COL5A2, THBS1, COL5A1, SPARC, and COL4A1) were identified as core genes based on the ranking of scores (Table 6). Furthermore, we also constructed a PPI network of hub genes with 10 nodes and 44 edges (Figure 5A) and interaction network between hub genes and their related genes (Figure 5B). The clustering coefficient of the PPI network was 0.978, and PPI enrichment P-value was lower than 1.0e-16.

Table 6

Top 10 Hub Genes in the PPI Network Ranked by the Degree Method

Rank	Name	Score
1	FN1	42
2	COL3A1	32
3	COL1A2	31
4	BGN	27
5	THBS2	26
5	COL5A2	26
7	THBS1	24
7	COL5A1	24
9	SPARC	23
10	COL4A1	22

Figure 5

(A) The interaction network of hub genes and their related genes. Red represents selected genes, and greed represents interacting genes. (B) PPI network of the hub genes constructed via the STRING online database.

As can be seen from the forest map (Figure 6A) and survival curve (Figure 6B–K) we constructed, among the 10 hub genes, 7 genes (FN1, COL1A2, BGN, THBS2, COL5A1, SPARC, COL4A1) were significantly associated with the OS of GC patients (Figure 6). High expression levels of all hub genes except COL3A1, COL5A2, and THBS1, were significantly related to poor OS in GC patients.

Figure 6

(A) A forest map of the 10 hub genes in all stages. 7 hub genes (FN1, COL1A2, BGN, THBS2, COL5A1, SPARC, and COL4A1) were significantly associated with the OS of GC patients. (B–K) The Kaplan-Meier survival curves of 10 hub genes.

Subsequently, genetic alteration information of the 10 hub genes was determined via cBioPortal, as illustrated in Figure 7. The genes COL1A2 and COL4A1 were altered most often (14% and 12%, respectively). These alterations included amplifications, missense mutations, deep deletions, truncating mutations, inframe mutations, and splice-site mutations (Figure 7A). Gene amplification is the most common of all types of genetic alterations (Figure 7B).

Figure 7

(A) Genetic alteration information of the 10 hub genes. (B) Alterations in the 10 hub genes from datasets of GC in the TCGA database.

Additionally, we verified the relative expression levels of hub genes using gene expression profiles of GC and normal tissues in TCGA and GTEx. As shown in Figure 8A–J, except THBS1, the expression levels of other hub genes were significantly up-regulated in GC compared with normal tissues (P < 0.05). The expression of THBS1 in GC tissues was up-regulated, but the difference was not statistically significant (P > 0.05). We further analyzed the expression levels of hub genes in different stages of GC (Figure 9). There were significant variations in the expression levels of BGN [Pr (>F) < 0.001] (Figure 9A), COL1A2 [Pr (>F) = 0.004] (Figure 9B), COL3A1 [Pr (>F) = 0.001] (Figure 9C), COL5A1 [Pr (>F) = 0.012] (Figure 9E), COL5A2 [Pr (>F) = 0.049] (Figure 9F), SPARC [Pr (>F) = 0.015] (Figure 9H), THBS1 [Pr (>F) = 0.011] (Figure 9I), and THBS2 [Pr (>F) = 0.025] (Figure 9J) in GC patients from different stages. These results strongly suggest that the hub genes identified hub genes may play an important role in the GC progression.

Figure 8

Comparison of expression levels of 10 hub genes in gastric cancer and normal tissues using the GEPIA2 database. We use red for gastric cancer tissues and gray for normal tissues. (A–J) respectively shows the BGN, COL1A2 COL3A1, COL4AI, COL5A1, COL5A2, FN1, SPARC, THBS1 and THBS2 expression level. Except for THBS1, the expression levels of the other nine hub genes in gastric cancer are significantly higher than those in normal tissues.

Figure 9

A comparison of expression levels of the 10 hub genes in different stages of GC. (A–J) respectively shows the BGN, COL1A2 COL3A1, COL4AI, COL5A1, COL5A2, FN1, SPARC, THBS1 and THBS2 expression level. There were significant variations in the expression levels of BGN (A), COL1A2 (B), COL3A1 (C), COL5A1 (E), COL5A2 (F), SPARC (H), THBS1 (I), and THBS2 (J) in GC patients from different stages.

Real-Time PCR and Immunohistochemistry

As shown in Figure 10A–N, immunohistochemical staining revealed that the expression of core genes in GC tissue was significantly higher than that in normal tissue. Statistical analysis confirmed that the relative expression levels in tumor tissues of mRNA of these 7 hub genes were significantly higher than in normal tissues (Figure 10O), and the average optical density (AOD) value of immunohistochemical staining was also significantly higher in tumor tissues (Figure 10P). The primer sequences used are shown in Table 7.

Figure 10

The results of real-time PCR and IHC. (A–N) IHC images of the 7 core genes in GC tissues and normal tissues under a 400 X microscope. (O) The relative hub gene mRNA levels in normal and GC tissues. (P) The AOD of hub genes expression. **P < 0.01; ***P < 0.001.

Table 7

Primer Sequences Used in PCR Experiments

Gene	Sequences
FN1-F	GGAACAAGCATGTCTCTCTGC
FN1-R	TGCTTGGAGAAGCTGTGAGT
COL1A2-F	ACTGTAAGAAAGGGCCCAGC
COL1A2-R	AGCAAAGTTCCCACCGAGAC
BGN-F	TGACTGGCATCCCCAAAGAC
BGN-R	GAGTAGCGAAGCAGGTCCTC
THBS2-F	ATCACACGCATCCGTCTCTG
THBS2-R	TAAACAGCCATTTGGGCATGG
COL5A1-F	TAGTCAGTACAAGCGGGGGTC
COL5A1-R	GTCATCTGTACCACACCCACA
SPARC-F	GCTTCGGCATCAAGCAGAGATA
SPARC-R	CATTGGGGGAAACACGAAGG
COL4A1-F	GGACCCCCGGGAGAAATAGG
COL4A1-R	GATTTGAAAAAGCAATGGCACTCC
GAPDH-F	GTCAGCCGCATCTTCTTT
GAPDH-R	CGCCCAATACGACCAAAT

A total of 348 GC patients were enrolled in the project. RNA sequencing data was available for 343 tumor tissue samples and 30 non-tumor tissue samples from 348 patients. We selected 7 hub genes with significant differences in survival analysis for diagnostic value evaluation. The ROC (receiver operating characteristic) curve (Figure 11) showed that these 7 hub genes had a high diagnostic value for GC: FN1 (AUC = 0.724, CI = 0.616–0.831), COL1A2 (AUC = 0.874, CI = 0.831–0.917), BGN (AUC = 0.949, CI = 0.925–0.973), THBS2 (AUC = 0.866, CI = 0.816–0.916), COL5A1 (AUC = 0.808, CI = 0.744–0.873), SPARC (AUC = 0.879, CI = 0.835–0.924), and COL4A1 (AUC = 0.861, CI = 0.783–0.939).

Figure 11

ROC curves of the 7 hub genes in the diagnosis of GC. Different colored curves represent the different core genes.

miRNA–mRNA Network

As illustrated in Figure 12, the miRNA-hub gene interaction network consists of 10 hub genes and 293 miRNAs. The four interactive hub genes that most miRNAs would target were THBS1 (degree score = 95), THBS2 (degree score = 81), COL4A1 (degree score = 42), and COL5A1 (degree score = 31). These were followed by COL1A2 (degree score = 30), SPARC (degree score = 21), COL3A1 (degree score = 11), BGN (degree score = 11), FN1 (degree score = 8), and COL5A2 (degree score = 7). Furthermore, hsa-miR-29b-3p (degree score = 6), hsa-miR-140-3p (degree, score = 5), hsa-miR-29a-3p (degree, score = 5), and hsa-miR-29c-3p (degree, score = 5) were the top four interactive miRNAs that targeted the most hub genes (Figure 12). Additionally, the detailed relationship between miRNAs and hub genes is shown in Table 8.

Figure 12

The interaction network between hub genes and target miRNAs. Hub genes are presented in yellow circles, whereas target miRNAs are shown in green circles. The interaction between hub genes and their related miRNAs is shown in the form of arrows.

Table 8

The Respective miRNAs Targeting the 10 Hub Genes

Gene	mi-RNA
BGN	hsa-miR-6758-5p, hsa-miR-6856-5p, hsa-miR-3125, hsa-miR-6859-5p, hsa-miR-4476, hsa-miR-3916, hsa-miR-4533, hsa-miR-6876-5p, hsa-miR-6828-5p, hsa-miR-8485, hsa-miR-4311
COL1A2	hsa-miR-4271, hsa-miR-6090, hsa-miR-4747-5p, hsa-miR-5196-5p, hsa-miR-4653-3p, hsa-miR-4424, hsa-miR-4781-5p, hsa-miR-6780b-5p, hsa-miR-4725-3p, hsa-miR-7111-5p, hsa-miR-4525, hsa-miR-6870-5p, hsa-miR-5698, hsa-miR-4723-5p, hsa-miR-625-5p, hsa-miR-1275, hsa-miR-569, hsa-miR-5010-5p, hsa-miR-6803-5p, hsa-miR-6751-5p, hsa-miR-7109-5p, hsa-miR-4665-5p, hsa-miR-4533, hsa-miR-26b-5p, hsa-miR-29c-3p, hsa-miR-1343-5p, hsa-miR-939-5p, hsa-miR-29a-3p, hsa-let-7 g-5p, hsa-miR-25-3p
COL3A1	hsa-miR-767-5p, hsa-miR-335-5p, hsa-miR-29c-3p, hsa-let-7b-5p, hsa-miR-143-3p, hsa-miR-455-3p, hsa-miR-29a-3p, hsa-miR-29a-5p, hsa-let-7d-5p, hsa-miR-29b-1-5p, hsa-miR-29b-3p
COL4A1	hsa-miR-218-5p, hsa-miR-124-3p, hsa-miR-767-5p, hsa-miR-3185, hsa-miR-542-3p, hsa-miR-4504, hsa-miR-5682, hsa-miR-6808-3p, hsa-miR-495-3p, hsa-miR-5688, hsa-miR-6871-3p, hsa-miR-562, hsa-miR-4716-5p, hsa-miR-21-5p, hsa-miR-7-1-3p, hsa-miR-7-2-3p, hsa-miR-125a-5p, hsa-miR-335-3p, hsa-miR-4318, hsa-miR-4693-5p, hsa-miR-7159-5p, hsa-miR-505-3p, hsa-miR-421, hsa-miR-8070, hsa-miR-498, hsa-miR-1305, hsa-miR-892a, hsa-miR-4272, hsa-miR-4709-5p, hsa-miR-29c-3p, hsa-miR-758-3p, hsa-miR-193b-3p, hsa-miR-18a-3p, hsa-miR-637, hsa-miR-769-3p, hsa-miR-16-5p, hsa-miR-92a-3p, hsa-miR-99a-5p, hsa-miR-100-5p, hsa-miR-29a-3p, hsa-miR-140-3p, hsa-miR-29b-3p
COL5A1	hsa-miR-3689a-3p, hsa-miR-30b-3p, hsa-miR-1273 h-5p, hsa-miR-6778-5p, hsa-miR-1233-5p, hsa-miR-6799-5p, hsa-miR-6883-5p, hsa-miR-6785-5p, hsa-miR-4728-5p, hsa-miR-6760-5p, hsa-miR-149-3p, hsa-miR-6878-5p, hsa-miR-7106-5p, hsa-miR-5703, hsa-miR-4516, hsa-miR-4434, hsa-miR-3934-5p, hsa-miR-764, hsa-miR-6780a-5p, hsa-miR-6779-5p, hsa-miR-3689c, hsa-miR-3689b-3p, hsa-miR-543, hsa-miR-26b-5p, hsa-miR-143-3p, hsa-miR-145-5p, hsa-miR-125a-3p, hsa-miR-6825-5p, hsa-miR-744-5p, hsa-miR-29b-1-5p, hsa-miR-29b-3p
COL5A2	hsa-miR-767-5p, hsa-miR-21-5p, hsa-miR-222-3p, hsa-miR-29c-3p, hsa-miR-143-3p, hsa-miR-29a-3p, hsa-miR-29b-3p
FN1	hsa-miR-218-5p, hsa-miR-200b-3p, hsa-miR-200c-3p, hsa-miR-26b-5p, hsa-miR-1-3p, hsa-miR-615-3p, hsa-let-7 g-5p, hsa-miR-140-3p
SPARC	hsa-miR-432-5p, hsa-miR-767-5p, hsa-miR-10a-5p, hsa-miR-192-5p, hsa-miR-29c-3p, hsa-miR-297, hsa-miR-675-3p, hsa-miR-3149, hsa-miR-2681-3p, hsa-miR-646, hsa-miR-4710, hsa-miR-4792, hsa-miR-7855-5p, hsa-miR-3653-3p, hsa-miR-4693-3p, hsa-miR-591, hsa-miR-29a-3p, hsa-miR-29a-5p, hsa-miR-433-5p, hsa-miR-211-5p, hsa-miR-29b-3p
THBS1	hsa-miR-3163, hsa-miR-548c-3p, hsa-miR-132-3p, hsa-miR-182-5p, hsa-miR-27b-3p, hsa-miR-155-5p, hsa-miR-19a-3p, hsa-miR-17-5p, hsa-miR-20a-5p, hsa-let-7a-5p, hsa-miR-98-5p, hsa-miR-3065-5p, hsa-miR-380-3p, hsa-miR-5003-5p, hsa-miR-202-3p, hsa-miR-567, hsa-miR-5096, hsa-miR-4500, hsa-let-7i-5p, hsa-let-7e-5p, hsa-miR-379-3p, hsa-miR-411-3p, hsa-miR-4458, hsa-let-7c-5p, hsa-let-7f-5p, hsa-miR-1252-3p, hsa-miR-7-1-3p, hsa-miR-7-2-3p, hsa-miR-4495, hsa-miR-3529-3p, hsa-miR-30a-5p, hsa-miR-30a-3p, hsa-miR-1-3p, hsa-let-7b-5p, hsa-miR-297, hsa-miR-675-3p, hsa-miR-3149, hsa-miR-92a-3p, hsa-miR-410-3p, hsa-miR-648, hsa-miR-6733-3p, hsa-miR-190a-3p, hsa-miR-4528, hsa-miR-4302, hsa-miR-4276, hsa-miR-4418, hsa-miR-509-5p, hsa-miR-509-3-5p, hsa-miR-211-3p, hsa-miR-4756-5p, hsa-miR-4739, hsa-miR-1321, hsa-miR-125a-3p, hsa-miR-6754-5p, hsa-miR-4441, hsa-miR-4270, hsa-miR-3665, hsa-miR-1205, hsa-miR-1273 g-3p, hsa-miR-6756-5p, hsa-miR-4450, hsa-miR-766-5p, hsa-miR-657, hsa-miR-198, hsa-miR-92a-2-5p, hsa-miR-6766-5p, hsa-miR-19b-3p, hsa-miR-3609, hsa-miR-6857-5p, hsa-miR-491-5p, hsa-miR-5680, hsa-miR-4637, hsa-miR-8084, hsa-miR-3121-3p, hsa-miR-1827, hsa-miR-377-5p, hsa-miR-4537, hsa-miR-6499-3p, hsa-miR-3612, hsa-miR-6086, hsa-miR-2467-3p, hsa-miR-455-3p, hsa-miR-650, hsa-miR-6825-5p, hsa-miR-6516-5p, hsa-miR-4655-5p, hsa-miR-4257, hsa-miR-6777-5p, hsa-miR-1910-3p, hsa-miR-6511a-5p, hsa-miR-6889-5p, hsa-let-7d-5p, hsa-let-7 g-5p, hsa-miR-487b-3p, hsa-miR-320c
THBS2	hsa-miR-27b-3p, hsa-miR-4719, hsa-miR-6814-5p, hsa-miR-662, hsa-miR-6793-3p, hsa-miR-4281, hsa-miR-133a-3p, hsa-miR-365a-3p, hsa-miR-365b-3p, hsa-miR-105-5p, hsa-miR-7853-5p, hsa-miR-6839-5p, hsa-miR-3183, hsa-miR-4667-3p, hsa-miR-6756-3p, hsa-miR-1470, hsa-miR-361-3p, hsa-miR-3127-3p, hsa-miR-378j, hsa-miR-6769b-3p, hsa-miR-4469, hsa-miR-4723-3p, hsa-miR-3191-5p, hsa-miR-1236-3p, hsa-miR-6826-3p, hsa-miR-6824-3p, hsa-miR-6515-3p, hsa-miR-518c-5p, hsa-miR-6764-3p, hsa-miR-4782-5p, hsa-miR-5706, hsa-miR-6767-5p, hsa-miR-6798-5p, hsa-miR-510-5p, hsa-miR-670-3p, hsa-miR-1343-5p, hsa-miR-4771, hsa-miR-4518, hsa-miR-1266-5p, hsa-miR-6762-5p, hsa-miR-939-5p, hsa-miR-6845-5p, hsa-miR-3922-5p, hsa-miR-4326, hsa-miR-3193, hsa-miR-1228-3p, hsa-miR-5572, hsa-miR-4491, hsa-miR-129-2-3p, hsa-miR-129-1-3p, hsa-miR-938, hsa-miR-299-5p, hsa-miR-4657, hsa-miR-6812-5p, hsa-miR-6737-5p, hsa-miR-4536-5p, hsa-miR-193b-5p, hsa-miR-6855-5p, hsa-miR-3170, hsa-miR-6819-5p, hsa-miR-6895-3p, hsa-miR-6873-3p, hsa-miR-593-3p, hsa-miR-6844, hsa-miR-4753-3p, hsa-miR-4768-5p, hsa-miR-6818-3p, hsa-miR-6809-3p, hsa-miR-6833-3p, hsa-miR-1307-3p, hsa-miR-6741-3p, hsa-miR-6845-3p, hsa-miR-130b-5p, hsa-miR-6817-3p, hsa-miR-3615, hsa-miR-1908-5p, hsa-miR-744-5p, hsa-miR-6787-5p, hsa-miR-663a, hsa-miR-135b-5p, hsa-miR-29b-3p

Drug-Gene Interaction Analysis

A total of 14 drugs related to 8 genes were selected, and the relationships between them were visualized through Cytoscape (Figure 13). Based on the total score on the website, the top ten drugs with higher scores were selected. The results are shown in Table 9.

Figure 13

The 14 predicted drugs associated with the hub genes. Red represents hub genes and blue represents potential drugs.

Table 9

Top 10 Drugs by Score Ranking in DGIdb

Gene	Drug	Interaction Type and Directionality	Interaction Score
FN1	L19IL2	n/a	15.95
FN1	L19TNFA	n/a	15.95
FN1	L19SIP 131I	n/a	15.95
THBS1	METHYLCELLULOSE	n/a	7.97
THBS2	BEVACIZUMAB	n/a	3.87
COL3A1	COLLAGENASE CLOSTRIDIUM HISTOLYTICUM	n/a	2.28
COL1A2	COLLAGENASE CLOSTRIDIUM HISTOLYTICUM	n/a	2.28
THBS1	VESNARINONE	n/a	2.28
THBS1	ABT-510	n/a	1.99
FN1	OCRIPLASMIN	Cleavage (inhibitory)	1.29

Discussion

Gastric Cancer remains the fifth most frequently diagnosed neoplasm and the third most deadly tumor in the world.14 With the rapid developments in molecular biology, leading to novel treatments, the current survival rate of GC patients has improved. However, there are certain limitations to treatment, and the outcome of patients is not always promising. Therefore, it is important to discover novel biomarkers and better understand the underlying molecular mechanisms to enable earlier diagnosis, treatment, and a better prognosis for GC patients. In recent years, with the rapid developments in the field of bioinformatics, more and more microarray and sequencing data can be obtained from online platforms, which provides a feasible way to discover new tumor-related diagnostic and therapeutic biomarkers.

In this study, three datasets for GC and normal tissue were downloaded from the GEO online database, and the differentially expressed genes were analyzed. A total of 340 common DEGs were identified, including 94 up-regulated and 246 downregulated genes. Subsequently, functional enrichment analysis of the DEGs was performed using the DAVID online tool. The top 3 significant terms from the GO enrichment analysis showed that the up-regulated DEGs were involved in the extracellular matrix organization, collagen catabolic process, and cell adhesion in the BP category. It is interesting to note that about 90% percent of tumor-related deaths are due to metastasis.15 Collagen is the main component of the extracellular matrix and interstitial microenvironment. It provides a scaffold for the growth of tumor cells and induces their migration.16 Jang et al have reported that collagen synthesis increases significantly in GC tissue.17 Additionally, changes in cell-cell adhesion and cell-matrix adhesion promote the chances of metastasis of the tumor cell.18 Collagen components are qualitatively and quantitatively reorganized in the tumor microenvironment. Collagen width has been identified as a useful prognostic indicator for GC.19 For KEGG pathway enrichment analysis, the up-regulated DEGs were enriched in ECM-receptor interaction, focal adhesion, and protein digestion and absorption. The structural integrity of ECM is crucial to maintain the normal and healthy activity of cells.20 Several studies have identified that the ECM plays an important role in various biological processes associated with tumor development, such as escaping apoptosis, cell growth regulation, tumor angiogenesis, invasion of cells, and metastasis.21–24 Thus, studying the remodeling of ECM that occurs in GC patients could lead to the discovery of potential therapeutic targets.

The degree of connectivity of a gene in a PPI network reflects its association with GC. In this study, a PPI network of 340 DEGs was constructed that included 321 nodes and 664 edges, and 10 hub genes were identified based on the highest degrees of connectivity. Subsequently, survival analysis of these hub genes was done, and it demonstrated that 7 hub genes (FN1, COL1A2, BGN, THBS2, COL5A1, SPARC, COL4A1) were significantly associated with the OS of GC patients. Then, the expression levels of these hub genes in tumor and normal tissues were compared through the GEPIA2 database. The results were consistent with those from the GEO database. PCR and IHC also confirmed that, compared with normal tissues, the expression levels of all hub genes were significantly higher in GC tissue. It is reasonable to hypothesize that these hub genes play an important role in the occurrence and development of GC. The prognostic value of hub genes was also evaluated based on data from TCGA. The ROC curve showed that 7 hub genes had a high diagnostic value for GC. Although a large number of studies have already shown that these hub genes are associated with the diagnosis and prognosis of GC, their underlying molecular mechanisms have not yet been fully clarified.

The SPARC gene encodes a cysteine-rich acidic matrix-associated protein which is required for collagen formation and is also involved in extracellular matrix synthesis and promotion of changes in cell shape.25 Previous studies have shown that the SPARC protein is associated with a variety of cancers, including GC;26 It may affect GC metastasis by regulating changes in the tumor microenvironment.27 These findings support the accuracy of our bioinformatics analyses to a certain extent.

The FN1 gene encodes fibronectin, a glycoprotein found as either a dimer or a polymer on the cell surface and in the extracellular matrix. Studies have shown that FN1 is associated with a variety of tumor types and that its expression is higher in tumor tissues.28–30 The FN1 protein probably interacts with vascular endothelial growth factor A and plays an important role in the development of non-small cell lung cancer (NSCLC). Hence, it can be used as a diagnostic or therapeutic target for NSCLC patients.31 However, no relevant studies on FN1 have been done in GC patients to date. Our study shows that the downregulation of FN1 is associated with better overall survival in GC patients. Furthermore, we identified 8 potential miRNAs (hsa-miR-218-5p, hsa-miR-200b-3p, hsa-miR-200c-3p, hsa-miR-26b-5p, hsa-miR-1-3p, hsa-miR-615-3p, hsa-let-7 g-5p, and hsa-miR-140-3p) that might be targets for FN1.

The THBS1 gene encodes thromboplastin 1, which is involved in angiogenesis and tumor progression, and its increased expression significantly correlates with tumor differentiation.32 Elevated expression of THBS1 protein has been associated with tumor growth and lymph node metastasis in GC patients.33 Interestingly, in this study, Kaplan-Meier’s plots showed that THBS1 overexpression was not an adverse prognostic factor for overall survival in patients with GC. Therefore, further studies on the effect of THBS1 on GC are needed.

Thromboplastin 2 (THBS2) is a member of the thromboplastin family of proteins that mediates the cell-cell and cell-matrix interactions. Overexpression of THBS2 has been associated with poor survival in patients with lung adenocarcinoma.34 THBS2 expression was also significantly correlate with the TNM stage, AJCC stage, and clinical outcome (P < 0.05). It may also be an important prognostic indicator of colorectal cancer.35 Previous studies have found that the expressions of TGF-β and THBS2 in GC cells are in a time-dependent relationship with the survival of patients,36 suggesting that THBS2 may affect PI3K by up-regulating the TGF-β signaling pathway.

BGN is a class I gene of the small-leucine-rich-proteoglycan (SLRP) family and encodes a leucine repeat-rich protein with a molecular weight of 42 kDa. BGN has been reported to be up-regulated in several cancer types and associated with tumor progression and worsening patient prognosis,37–39 which is consistent with our findings. As a potential new prognostic biomarker for GC, BGN can be easily translated into the clinical setting and provides an opportunity to develop new therapies.

The COL family of genes is involved in the interactions of the ECM-receptor and local adhesion pathway, which play an important role in the invasion and metastasis of tumors. Collagen, laminin, and fibronectin are major components of the extracellular matrix, which provides the matrix for cancer cell growth, survival, and migration.40 Studies have shown that the COL1A2 mRNA levels are related to the histological types and lymph node status of GC, and may be involved in the invasion, metastasis, and carcinogenesis of GC.41,42 All of these previous findings support the accuracy of our results. In addition, the expression of COL5A2 is also closely related to the occurrence and development of colon cancer.43 Overexpression of COL3A1 has been confirmed in a variety of cancers, such as bladder cancer,44,45 but the effect of its expression levels in GC is not fully understood. COL5A1 has been reported to be closely associated with breast cancer,46 lung adenocarcinoma,47 oral squamous cell carcinoma,48 and GC.49 High expression of COL4A1 was found to be associated with an advanced tumor stage and poor overall survival in patients with hepatocellular carcinoma.50 Furthermore, in breast cancer cells, COL4A1 knockout resulted in reduced cell viability and cell cycle arrest.51 In this study, we found that COL1A2, COL3A1, COL4A1, COL5A1, and COL5A2 were significantly overexpressed in GC tissue, which is consistent with the TCGA results, indicating that the COL family of genes could make good biomarker panels for GC diagnosis.

With the help of the DGIdb database, a total of 14 drugs with potential therapeutic effects on GC were identified. The top 10 drugs among these, based on the website’s total score, were L19IL2, L19TNFA, L19SIP 131I, methylcellulose, bevacizumab, collagenase clostridium histolyticum, vesnarinone, ABT-510, and ocriplasmin. These drugs were mainly related to genes FN1, COL3A1, COL1A2, THBS2, COL5A2, THBS1, COL5A1, and COL4A1. The ability of the L19 antibody to preferentially target tumor sites has been widely demonstrated using nuclear medicine procedures with the preparation of radiologically labeled antibodies. L19IL2 is a recombinant fusion protein, containing the L19 antibody in diabody format fused to human IL2.52 In a study of the treatment of advanced metastatic melanoma, it was found that the use of L19IL2 in combination with dacarbazine (DTIC) significantly improved the overall response rates and median progression-free survival, as compared to the use of DTIC alone.53 The use of radioimmunotherapy with 131I-L19SIP has been reported to induce a sustained partial response in patients with relapsed Hodgkin’s lymphoma.54 L19-TNFα/L19-IL2 combination therapy could eradicate tumors in myeloma mice through TNFα-induced tumor necrosis and L19-TNFα/L19-IL2-mediated local cellular immune response.55 Surfactant-free solubilization of docetaxel is an efficient method for insoluble drug delivery to enhance the antitumor effect.56 Bevacizumab is expected to remain a key agent in cancer therapy, both due to its established efficacy in approved indications and its promise as a partner in novel targeted combination treatments.57 Combined bevacizumab can significantly reduce the risk of progression of advanced tumors and improve the overall survival rate in cases of metastatic colorectal cancer,58 non-small cell lung cancer,59 and renal cell carcinoma.60 Vesnarinone is a unique anti-proliferating, differentiation-inducing, and apoptosis-inducing drug against several human malignancies, including leukemia and several solid tumors.61 ABT-510 is a small mimetic peptide of a natural anti-angiogenic compound called thromboreactive protein-1. It has been shown to significantly reduce tumor growth and burden in preclinical mouse models and in naturally occurring dog tumors.62 Additionally, ABT-510 induced tumor cell apoptosis and inhibited ovarian tumor growth in an orthotopic, syngeneic model of epithelial ovarian cancer.62 It is interesting to note that all the above-mentioned drugs are gene-related and obtained through database predictions. They provide us with new potential strategies for the treatment of cancers. However, to determine whether they can be used in GC treatment, we need more experiments and clinical verification.

Conclusions and Limitations

In this study, we have successfully screened 10 hub genes associated with GC using the GEO database and bioinformatics analyses and have established that 7 of these genes were associated with poor cancer prognosis. We have also studied the mRNA and protein expression levels through real-time PCR and IHC experiments. We demonstrated that searching for drugs related to these hub genes in the DGIdb has provided new insights for GC treatment. However, this study has some limitations. Due to the lack of further experimental studies and validation, we cannot explain exactly how the hub gene-miRNAs networks and therapeutic drugs play a role in the diagnosis and treatment of GC. Due to the short follow-up time and a small number of patients (n = 60), clinicopathological data could not be analyzed conclusively. This issue will be addressed in future research.

In summary, identifying the hub genes, miRNA–mRNA network, and potential candidate drugs associated with GC provides new insights into the molecular mechanisms and treatment of GC.