Identification of circRNA-miRNA-Immune-Related mRNA Regulatory Network in Gastric Cancer.

The pathogenesis of gastric cancer (GC) is still not fully understood. We aimed to find the potential regulatory network for ceRNA (circRNA-miRNA-immune-related mRNA) to uncover the pathological molecular mechanisms of GC. The expression profiles of circRNA, miRNA, and mRNA in gastric tissue from GC patients were downloaded from the Gene Expression Omnibus (GEO) datasets. Differentially expressed circRNAs, miRNAs, and immune-related mRNAs were filtered, followed by the construction of the ceRNA (circRNA-miRNA-immune-related mRNA) network. Functional annotation and protein-protein interaction (PPI) analysis of immune-related mRNAs in the network were performed. Expression validation of circRNAs and immune-related mRNAs was performed in the new GEO and TCGA datasets and in-vitro experiment. A total of 144 differentially expressed circRNAs, 216 differentially expressed miRNAs, and 2,392 differentially expressed mRNAs were identified in GC. Some regulatory pairs of circRNA-miRNA-immune-related mRNA were obtained, including hsa_circ_0050102-hsa-miR-4537-NRAS-Tgd cells, hsa_circ_0001013-hsa-miR-485-3p-MAP2K1-Tgd cells, hsa_circ_0003763-hsa-miR-145-5p-FGF10-StromaScore, hsa_circ_0001789-hsa-miR-1269b-MET-adipocytes, hsa_circ_0040573-hsa-miR-3686-RAC1-Tgd cells, and hsa_circ_0006089-hsa-miR-5584-3p-LYN-neurons. Interestingly, FGF10, MET, NRAS, RAC1, MAP2K1, and LYN had potential diagnostic value for GC patients. In the KEGG analysis, some signaling pathways were identified, such as Rap1 and Ras signaling pathways (involved NRAS and FGF10), Fc gamma R-mediated phagocytosis and cAMP signaling pathway (involved RAC1), proteoglycans in cancer (involved MET), T-cell receptor signaling pathway (involved MAP2K1), and chemokine signaling pathway (involved LYN). The expression validation of hsa_circ_0003763, hsa_circ_0004928, hsa_circ_0040573, FGF10, MET, NRAS, RAC1, MAP2K1, and LYN was consistent with the integrated analysis. In conclusion, the identified ceRNA (circRNA-miRNA-immune-related mRNA) regulatory network may be associated with the development of GC.

Introduction

Gastric cancer (GC) is one of the most serious malignant tumors (1). According to the GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries, the incidence of GC is 7.7% (2). In addition, with 2.26 million new cases estimated in 2020, GC (0.77 million) has become the most commonly diagnosed cancer worldwide, ranking only second to female breast cancer (3). Recurrence is the main cause of GC-related death. Although mortality is steadily decreasing, GC still leads to a poor diagnosis and prognosis for patients (4). In addition, the 5-year survival rate of patients is still very low in serious GC patients. It has been shown that obesity, active tobacco smoking, high meat and salt intake, low vegetable/fruit intake, Helicobacter pylori infection, and gut microbiota have been shown to be associated with an increased risk of GC (5–8). In addition, epigenetic alterations are associated with the processes of gastric carcinogenesis and metastasis (9). Clinically, surgery is the only curative treatment; however, some patients have inoperable disease at diagnosis (10). Hence, there is a need to elucidate the potential molecular mechanisms in the development of GC and to look for new molecular markers and therapeutic targets.

circRNAs usually result from splicing/back-splicing events via exon or intron circularization (11). circRNA can act as a molecular sponge of miRNA to regulate mRNA expression. It was reported that sponge circRNAs are part of a complex RNA-binding protein–circRNA–miRNA–mRNA interaction network and are involved in the establishment, chemoresistance, and progression of GC (12). hsa_circ_LARP4 can inhibit cell invasion of GC by sponging hsa-miR-424-5p and regulating large tumor suppressor kinase 1 (LATS1) expression (13). hsa_circ_NRIP1 can act as an hsa-miR-149-5p sponge to promote GC progression via the AKT serine/threonine kinase 1 (AKT1)/mechanistic target of rapamycin kinase (mTOR) pathway (14). hsa_circ_NHSL1 can promote GC progression through the hsa-miR-1306-3p/SIX homeobox 1 (SIX1)/vimentin axis (15). hsa_circ_CACTIN can promote GC progression by sponging hsa-miR-331-3p and regulating transforming growth factor beta receptor 1 (TGFBR1) expression (16). hsa_circ_0026359 can enhance cisplatin resistance in GC via targeting the hsa-miR-1200/DNA polymerase delta 4, accessory subunit (POLD4) pathway (17). In addition, it is believed that the pathogenies and progression of GC are influenced by the cross-talk between tumor cells and the host immune system (18–20). In view of this, we tried to find the potential differentially expressed circRNAs, miRNAs, and immune-related mRNAs in GC.

Materials and Methods

Data Retrieval and Analysis

In this study, the expression profiles of circRNA, miRNA, and mRNA were downloaded from the Gene Expression Omnibus (GEO) datasets by searching keywords [“gastric cancer” (All Fields) AND “Homo sapiens”(porgn) AND “gse”(Filter)]. The following datasets were selected: 1) dataset must be genome-wide transcriptome data of mRNA/miRNA/circRNA; 2) data were obtained from tumor tissues of the GC group and paracancer control group; and 3) both standardized and raw datasets were considered. Finally, two circRNA expression datasets (GSE83521 and GSE89143), two miRNA expression datasets (GSE93415 and GSE158315), and two mRNA expression datasets (GSE66229 and GSE65801) were selected ( ).

Table 1

Detailed information of selected datasets of circRNA, miRNA, and mRNA expression datasets.

	GEO accession	Author	Platform	Samples (N:GC)	Year	Tissue
circRNA	GSE83521	Yan Zhang	GPL19978 Agilent-069978 Arraystar Human CircRNA microarray V1	6:6	2017	Gastric tissue
	GSE89143	Junming Guo	GPL19978 Agilent-069978 Arraystar Human CircRNA microarray V1	3:3	2017	Gastric tissue
miRNA	GSE93415	Marek Sierżęga	GPL19071 Exiqon miRCURY LNA microRNA array; 7th generation REV - hsa, mmu, and rno; batch 208520-22; lot 35101-35101 (miRBase 19.0)	20:20	2017	Gastric tissue
	GSE158315	Yuping Wang	GPL18058 Exiqon miRCURY LNA microRNA array, 7th generation (miRBase v18, condensed Probe_ID version)	5:5	2021	Gastric tissue
mRNA	GSE66229	Michael Nebozhyn	GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array	100:300	2015	Gastric tissue
	GSE65801	Hao Li	GPL14550 Agilent-028004 SurePrint G3 Human GE 8 × 60K Microarray (Probe Name Version)	32:32	2015	Gastric tissue

N, paracancer control group; GC, gastric cancer group.

Screening of Differentially Expressed circRNAs, miRNAs, and Immune-Related mRNAs

The probe and ID of circRNA/miRNA/mRNA were mapped one by one. After scale standardization, the dataset was merged and batch effects were removed via the ComBat function of sva package (R-4.0.5). The metaMA and limma packages were used to identify circRNAs/miRNAs/mRNAs. The default parameter of the Pvalcombination command of the metaMA package was used to make the difference. p-values and effect sizes (ES, the effectSize obtained from the metaMA package) from data were calculated either from classical or moderated t-tests. These p-values were combined by the inverse normal method. The Benjamini–Hochberg threshold was used to calculate the false discovery rate (FDR). |Combined.ES| >1 and FDR <0.05 were the screening criteria for circRNAs/miRNAs/mRNAs. In addition, immune-related mRNAs were downloaded from the ImmPort database (https://www.immport.org/shared/home). Those immune-related mRNAs were obtained by intersection of differentially expressed mRNAs and immune-related mRNAs in the ImmPort database. xCell (21) was used to calculate the distribution of immune cells in each sample based on the ssGSEA method. xCell score was organized into immune cell infiltration matrix to calculate the types of immune cells that differed between the GC group and the normal control group. The Pearson correlation coefficient method was used to calculate the correlation between mRNAs and differential immune cells.

Construction of the ceRNA (circRNA–miRNA–Immune-Related mRNA) Regulatory Network

The TargetScan (http://www.targetscan.org/vert_71/) software was utilized to predict targeted relationship between differentially expressed circRNAs and differentially expressed miRNAs. In addition, the miRWalk software (http://mirwalk.umm.uni-heidelberg.de/interactions/) was applied to predict targeted relationship between miRNAs and mRNAs. The relationship pairs of miRNA–mRNA verified in at least one database (TargetScan, miRDB, and MiRTarBase) were selected. The overlapping mRNAs were obtained between predicted mRNAs (in the miRNA–mRNA relationship pairs) and immune-related differential expression mRNAs. The ceRNA (circRNA–miRNA–immune-related mRNA) regulatory network was constructed by fusing with the circRNA–miRNA relationship pairs and miRNA–immune-related mRNA relationship pairs.

Functional Annotation of Immune-Related mRNAs in the ceRNA Regulatory Network

To investigate the function of immune-related mRNAs in the regulatory ceRNA network, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed by using the DAVID database (https://david.ncifcrf.gov/tools.jsp). FDR <0.05 was considered as statistical significance.

Protein–Protein Interaction Network of Immune-Related mRNAs in the ceRNA Regulatory Network

To further explore the interaction between immune-related mRNAs in the ceRNA regulatory network, protein–protein interaction (PPI) was performed by using the STRING database. The results were imported using the Cytoscape software (http://www.cytoscape.org/). The CytoHubba plug-in was used to filter core immune-related mRNAs by intersecting the first 10 mRNAs of each algorithm (degree, MNC, MCC, and EPC). In addition, the ROC analysis was carried out to assess the diagnostic value of core immune-related mRNAs in the PPI network.

Expression Validation of Differentially Expressed circRNAs and Immune-Related mRNAs

In order to further validate the expression of identified circRNAs and immune-related mRNAs, electronic validation was performed. The GSE93541 dataset (involving tumor tissues from three cases and three normal controls) and the GSE141977 dataset (involving plasma from three cases and three normal controls) were used for expression validation of identified circRNAs. Additionally, the Cancer Genome Atlas (TCGA) dataset (involved tumor tissues from 375 cases and 32 normal controls) was used to validate the expression of identified mRNAs. The expression result of these circRNAs and mRNAs was shown by the box plots. Statistical significance was ascribed to p-value <0.05.

In addition, qRT-PCR was used to further validate the expression of identified circRNAs and immune-related mRNAs. Five GC patients and five normal controls were enrolled. The blood sample of these individuals was collected. All participating individuals provided informed consent with the approval of the Ethics Committee of Zhejiang Hospital (2021-164K). Total RNA of the blood sample was extracted and DNA was synthesized using FastQuant cDNA first-strand synthesis kit. qRT-PCR was performed in the SuperReal PreMix Plus (SYBR Green). Relative circRNA/mRNA expression was analyzed by the 2−△△CT method. 2−△△CT >1 and 2−△△CT <1 represented upregulation and downregulation, respectively.

Results

Differentially Expressed circRNAs, miRNAs, and Immune-Related mRNAs

A total of 144 differentially expressed (55 upregulated and 89 downregulated) circRNAs, 216 differentially expressed (105 upregulated and 111 downregulated) miRNAs, and 2,392 differentially expressed (1,329 upregulated and 1,063 downregulated) mRNAs were identified in GC. The top 10 differentially expressed circRNAs were identified, such as hsa_circ_0045602 (upregulation), hsa_circ_0006089 (upregulation), hsa_circ_0001789 (upregulation), hsa_circ_0018004 (downregulation), and hsa_circ_0003763 (downregulation) ( ). hsa-miR-4537 was one of top 10 downregulated miRNAs ( ). The top 10 differentially expressed mRNAs are listed in . The heat map of the top 100 differentially expressed circRNAs, miRNAs, and mRNAs is shown in , respectively. In addition, 1,793 immune-related mRNAs were obtained from the ImmPort database. It was noted that 159 immune-related mRNAs were also differentially expressed mRNAs in GC ( ).

Table 2

The top 10 upregulated and downregulated circRNAs in GC.

ID	Symbol	Combined.ES	p-value	FDR	Up/down
hsa_circRNA_102191	hsa_circ_0045602	3.118562462	1.07E−06	0.000830662	Up
hsa_circRNA_104947	hsa_circ_0007613	3.076836769	1.43E−06	0.000830662	Up
hsa_circRNA_101882	hsa_circ_0040573	3.859880212	1.82E−06	0.000830662	Up
hsa_circRNA_102592	hsa_circ_0052001	3.180131272	2.28E−06	0.000830662	Up
hsa_circRNA_101471	hsa_circ_0034398	2.865182594	3.78E−06	0.000836756	Up
hsa_circRNA_102614	hsa_circ_0006089	2.831849087	3.94E−06	0.000836756	Up
hsa_circRNA_104589	hsa_circ_0001789	3.306599675	4.20E−06	0.000836756	Up
hsa_circRNA_100641	hsa_circ_0019054	2.89673611	4.77E−06	0.00087084	Up
hsa_circRNA_101875	hsa_circ_0040481	2.642927097	9.51E−06	0.001488245	Up
hsa_circRNA_102777	hsa_circ_0055521	2.42548772	2.78E−05	0.003039928	Up
hsa_circRNA_100571	hsa_circ_0018004	−3.196809378	9.51E−07	0.000830662	Down
hsa_circRNA_104599	hsa_circ_0001793	−3.624031186	2.17E−06	0.000830662	Down
hsa_circRNA_400066	hsa_circ_0092330	−2.876804101	3.17E−06	0.000836756	Down
hsa_circRNA_101651	hsa_circ_0036941	−2.829679697	3.88E−06	0.000836756	Down
hsa_circRNA_001914	hsa_circ_0000902	−2.987023494	7.15E−06	0.001204509	Down
hsa_circRNA_101965	hsa_circ_0000740	−2.632340492	1.05E−05	0.001531677	Down
hsa_circRNA_103442	hsa_circ_0003763	−2.780691514	1.88E−05	0.00256875	Down
hsa_circRNA_400056	hsa_circ_0092297	−2.597537412	2.21E−05	0.00275069	Down
hsa_circRNA_100382	hsa_circ_0007277	−2.696712581	2.26E−05	0.00275069	Down
hsa_circRNA_104661	hsa_circ_0004366	−2.495771422	2.39E−05	0.00275519	Down

ES, effect size; FDR, false discovery rate.

Table 3

The top 10 upregulated and downregulated miRNAs in GC.

Symbol	Combined.ES	p-value	FDR	Up/down
hsa-miR-181a-5p	2.621880675	8.01E−13	6.84E−10	Up
hsa-miR-181b-5p	2.541945488	4.05E−12	1.15E−09	Up
hsa-miR-23a-3p	2.283816919	1.17E−10	1.11E−08	Up
hsa-miR-331-3p	2.121044968	1.04E−09	6.82E−08	Up
hsa-miR-320c	2.029389583	1.53E−09	8.73E−08	Up
hsa-miR-25-3p	1.94259155	4.78E−09	2.55E−07	Up
hsa-miR-320d	1.961475455	6.29E−09	2.98E−07	Up
hsa-miR-92a-3p	1.928060578	7.06E−09	3.02E−07	Up
hsa-miR-21-3p	1.858793042	1.58E−08	5.18E−07	Up
hsa-miR-320e	1.904827788	1.72E−08	5.44E−07	Up
hsa-miR-4279	−2.564714404	1.76E−12	7.51E−10	Down
hsa-miR-3124-3p	−2.474608778	5.60E−12	1.19E−09	Down
hsa-miR-4728-3p	−2.455455298	1.21E−11	2.07E−09	Down
hsa-miR-4635	−2.393092709	1.54E−11	2.19E−09	Down
hsa-miR-4537	−2.418295745	5.96E−11	7.27E−09	Down
hsa-miR-642b-5p	−2.298345262	8.63E−11	9.21E−09	Down
hsa-miR-5196-3p	−2.172908344	3.75E−10	3.20E−08	Down
hsa-miR-877-3p	−2.202854847	4.82E−10	3.51E−08	Down
hsa-miR-4290	−2.27123672	4.93E−10	3.51E−08	Down
hsa-miR-4268	−2.219462893	1.40E−09	8.54E−08	Down

ES, effect size; FDR, false discovery rate.

Table 4

The top 10 upregulated and downregulated mRNAs in GC.

ID	Symbol	Combined.ES	p-value	FDR	Up/down
9997	SCO2	1.400471	<0.001	<0.001	Up
9991	ROD1	1.325606	<0.001	<0.001	Up
9966	TNFSF15	1.160469	<0.001	<0.001	Up
995	CDC25C	1.639835	<0.001	<0.001	Up
994	CDC25B	1.62622	<0.001	<0.001	Up
993	CDC25A	1.407589	<0.001	<0.001	Up
9928	KIF14	2.287231	<0.001	<0.001	Up
9926	LPGAT1	1.19746	<0.001	<0.001	Up
9918	NCAPD2	1.645622	<0.001	<0.001	Up
991	CDC20	1.774075	<0.001	<0.001	Up
9992	KCNE2	−1.72329	<0.001	<0.001	Down
9934	P2RY14	−1.90515	<0.001	<0.001	Down
9905	SGSM2	−1.00543	<0.001	<0.001	Down
9892	SNAP91	−1.27282	<0.001	<0.001	Down
9886	RHOBTB1	−1.07693	<0.001	<0.001	Down
9874	TLK1	−1.11589	<0.001	<0.001	Down
9867	PJA2	−1.7726	<0.001	<0.001	Down
9832	JAKMIP2	−1.09529	<0.001	<0.001	Down
9829	DNAJC6	−1.03815	<0.001	<0.001	Down
9783	RIMS3	−1.64884	<0.001	<0.001	Down

ES, effect size; FDR, false discovery rate.

Figure 1

The heat map of the top 100 differentially expressed circRNAs (A), miRNAs (B), and mRNAs (C) in gastric cancer (GC).

Figure 2

Venn diagram of differentially expressed mRNAs in GC and immune-related mRNAs in the ImmPort database. dEmRNA, differentially expressed mRNA.

ceRNA (circRNA–miRNA–Immune-Related mRNA) Regulatory Network

A total of 2,037 negatively regulated targeting relationship pairs of circRNA–miRNA (involving 143 circRNAs and 198 miRNAs) and 142 negatively regulated targeting relationship pairs of miRNA–immune-related mRNA (involving 97 miRNAs and 58 immune-related mRNAs) were identified via the TargetScan software and miRWalk, respectively. The ceRNA (circRNA–miRNA–immune-related mRNA) regulatory network (involving 137 circRNAs, 96 miRNAs, and 58 immune-related mRNAs) was obtained by fusing with the circRNA–miRNA relationship pairs and miRNA–immune-related mRNA relationship pairs ( ). Some ceRNA relationship pairs were identified, such as hsa_circ_0045602–hsa-miR-4538–WNT5A and hsa_circ_0018004–hsa-miR-199a-5p–KL.

Figure 3

The ceRNA (circRNA–miRNA–immune-related mRNA) regulatory network in GC. Circle, V shape, and rectangle represent differentially expressed circRNA, miRNA, and mRNA, respectively. Red and green represent upregulation and downregulation, respectively. The black border represents the top 10 up/downregulated circRNAs/miRNAs/mRNAs.

Enrichment Analysis of Immune-Related mRNAs in the ceRNA Regulatory Network

Fifty-eight immune-related mRNAs in the regulatory ceRNA network were used for functional enrichment analysis. GO analysis showed that signal transduction, extracellular space, and growth factor activity were the most significantly enriched biological process, cytological component, and molecular function, respectively ( ). According to the KEGG analysis, Rap1 and Ras signaling pathways [involved NRAS proto-oncogene, GTPase (NRAS), and fibroblast growth factor 10 (FGF10)], Fc gamma R-mediated phagocytosis and cAMP signaling pathway [involved Rac family small GTPase 1 (RAC1)], proteoglycans in cancer [involved MET proto-oncogene, receptor tyrosine kinase (MET)], T-cell receptor signaling pathway [involved mitogen-activated protein kinase kinase 1 (MAP2K1)], and chemokine signaling pathway [involved LYN proto-oncogene, Src family tyrosine kinase (LYN)] were significantly enriched signaling pathways ( and ).

Figure 4

GO (A) and KEGG (B) enrichment analyses of immune-related mRNAs in the regulatory ceRNA network in GC. BP, biological process; CC, cytological component; MF, molecular function.

Table 5

KEGG enrichment analysis of immune-related mRNAs in the ceRNA regulatory network in GC.

Term	Count	p-value	FDR	mRNAs
hsa05200:Pathways in cancer	16	8.65E−09	9.34E−07	MAP2K1, PRKCB, F2R, WNT5A, TGFBR1, PGF, NRAS, FGF14, TRAF3, CDK4, AGTR1, RAC1, FGF13, BID, MET, FGF10
hsa04060:Cytokine–cytokine receptor interaction	11	2.17E−06	1.17E−04	GHR, CXCL9, TNFRSF12A, TNFSF15, TNFRSF9, LIF, TNFSF11, TNFRSF10B, LIFR, IL1RAP, TGFBR1
hsa05218:Melanoma	7	5.05E−06	1.47E−04	MAP2K1, NRAS, FGF14, CDK4, FGF13, MET, FGF10
hsa04015:Rap1 signaling pathway	10	5.43E−06	1.47E−04	MAP2K1, NRAS, FGF14, PRKCB, F2R, RAC1, FGF13, MET, PGF, FGF10
hsa04014:Ras signaling pathway	10	9.89E−06	2.14E−04	MAP2K1, NRAS, FGF14, PRKCB, RAC1, FGF13, PAK3, MET, PGF, FGF10
hsa04810:Regulation of actin cytoskeleton	9	4.54E−05	7.33E−04	MAP2K1, NRAS, FGF14, F2R, RAC1, FGF13, PAK3, VAV2, FGF10
hsa04151:PI3K–Akt signaling pathway	11	4.75E−05	7.33E−04	GHR, MAP2K1, NRAS, FGF14, CDK4, F2R, RAC1, FGF13, MET, PGF, FGF10
hsa04664:Fc epsilon RI signaling pathway	6	6.55E−05	8.85E−04	LYN, MAP2K1, NRAS, PRKCB, RAC1, VAV2
hsa04650:Natural killer cell mediated cytotoxicity	7	1.12E−04	0.001343	MAP2K1, NRAS, PRKCB, TNFRSF10B, RAC1, BID, VAV2
hsa05211:Renal cell carcinoma	5	7.86E−04	0.008494	MAP2K1, NRAS, RAC1, PAK3, MET
hsa04662:B-cell receptor signaling pathway	5	9.30E−04	0.009003	LYN, MAP2K1, NRAS, RAC1, VAV2
hsa04010:MAPK signaling pathway	8	0.001	0.009003	MAP2K1, NRAS, FGF14, PRKCB, RAC1, FGF13, TGFBR1, FGF10
hsa04612:Antigen processing and presentation	5	0.001337	0.011104	PSME3, TAP2, TAP1, CALR, CTSB
hsa05205:Proteoglycans in cancer	7	0.001573	0.012138	MAP2K1, NRAS, PRKCB, WNT5A, PLAUR, RAC1, MET
hsa04510:Focal adhesion	7	0.00183	0.013074	MAP2K1, PRKCB, RAC1, PAK3, MET, PGF, VAV2
hsa04666:Fc gamma R-mediated phagocytosis	5	0.001937	0.013074	LYN, MAP2K1, PRKCB, RAC1, VAV2
hsa04660:T-cell receptor signaling pathway	5	0.003655	0.023221	MAP2K1, NRAS, CDK4, PAK3, VAV2
hsa05223:Non-small cell lung cancer	4	0.005345	0.032068	MAP2K1, NRAS, PRKCB, CDK4
hsa04062:Chemokine signaling pathway	6	0.006326	0.034272	LYN, MAP2K1, CXCL9, NRAS, RAC1, VAV2
hsa04730:Long-term depression	4	0.006482	0.034272	LYN, MAP2K1, NRAS, PRKCB
hsa04370:VEGF signaling pathway	4	0.006787	0.034272	MAP2K1, NRAS, PRKCB, RAC1
hsa04071:Sphingolipid signaling pathway	5	0.006981	0.034272	MAP2K1, NRAS, PRKCB, RAC1, BID
hsa05214:Glioma	4	0.008093	0.035418	MAP2K1, NRAS, PRKCB, CDK4
hsa05212:Pancreatic cancer	4	0.008093	0.035418	MAP2K1, CDK4, RAC1, TGFBR1
hsa04024:cAMP signaling pathway	6	0.008199	0.035418	MAP2K1, FSHB, F2R, RAC1, ADRB2, VAV2
hsa05220:Chronic myeloid leukemia	4	0.010713	0.044498	MAP2K1, NRAS, CDK4, TGFBR1
hsa04550:Signaling pathways regulating pluripotency of stem cells	5	0.01188	0.047521	MAP2K1, NRAS, WNT5A, LIF, LIFR

FDR, false discovery rate.

PPI Network of Immune-Related mRNAs in the ceRNA Regulatory Network

In order to further explore the interaction between immune-related mRNAs in the regulatory ceRNA network, the PPI network was established ( ). Six core immune-related mRNAs were identified by four algorithms, namely, FGF10, MET, NRAS, RAC1, MAP2K1, and LYN ( ). It is worth mentioning that FGF10 (AUC = 0.764), MET (AUC = 0.829), NRAS (AUC = 0.840), RAC1 (AUC = 0.797), MAP2K1 (AUC = 0.813), and LYN (AUC = 0.817) had a potential diagnostic value for GC ( ). It was indicated that these mRNAs had a diagnostic potential for GC. In addition, the ceRNA subnetwork (involving 59 circRNAs and 10 miRNAs) based on the above six mRNAs was constructed ( ). The ceRNA subnetwork contains 75 nodes and 124 edges. Some ceRNA relationship pairs were identified, such as hsa_circ_0050102–hsa-miR-4537–NRAS, hsa_circ_0001013–hsa-miR-485-3p–MAP2K1, hsa_circ_0003763–hsa-miR-145-5p–FGF10, hsa_circ_0001789–hsa-miR-1269b–MET, hsa_circ_0040573–hsa-miR-3686–RAC1, and hsa_circ_0006089–hsa-miR-5584-3p–LYN. In addition, xCell was used to calculate the distribution of immune cells in each sample based on the ssGSEA method. The types of immune cells were different in the GC group from those in the normal control group. For example, Tgd cells and adipocytes were the most upregulated and downregulated types of immune cells in GC ( ). The Pearson correlation coefficient method was used to calculate the correlation between above six mRNAs of differential immune cells ( ). The result showed that NRAS, RAC1, and MAP2K1 were positively correlated with Tgd cells. FGF10 was positively correlated with StromaScore. LYN and MET were negatively correlated with neurons and adipocytes, respectively.

Figure 5

The PPI network of immune-related mRNAs in the regulatory ceRNA network in GC.

Figure 6

Identification of the six core immune-related mRNAs in the PPI network via four algorithms in GC.

Figure 7

ROC analysis of six core immune-related mRNAs in the PPI network in GC.

Figure 8

The ceRNA subnetwork based on six core immune-related mRNAs in the PPI network in GC. Circle, V shape, and rectangle represent differentially expressed circRNA, miRNA, and mRNA, respectively. Red and green represent upregulation and downregulation, respectively. The black border represents the top 10 up/downregulated circRNAs/miRNAs/mRNAs.

Figure 9

Distribution of immune cell types in GC. *p < 0.05; **p < 0.01; ***p < 0.001.

Figure 10

The correlation analysis between six core immune-related mRNAs and immune cells in GC. *p < 0.05; **p < 0.01; ***p < 0.001.

To further validate the expression of identified circRNAs and immune-related mRNAs, electronic validation was performed. The GSE93541 dataset ( ) and the GSE141977 dataset ( ) were used for expression validation of hsa_circ_0003763, hsa_circ_0004928, and hsa_circ_0040573. The result showed that hsa_circ_0003763 was downregulated and hsa_circ_0004928 and hsa_circ_0040573 were upregulated in GC. Although there was no significant difference (may be caused by the small sample size), the expression trend was consistent with the integrated analysis results. In addition, the TCGA dataset was used to validate the expression of FGF10, MET, NRAS, RAC1, MAP2K1, and LYN ( ). The result showed that FGF10 was significantly downregulated. MET, NRAS, RAC1, MAP2K1, and LYN were remarkably upregulated in GC, which is in line with the integrated analysis results.

Figure 11

Expression validation of hsa_circ_0003763, hsa_circ_0004928, and hsa_circ_0040573 in the GSE93541 dataset (A) and GSE141977 dataset (B). Ns, no significant difference.

Figure 12

Expression validation of FGF10, MET, NRAS, RAC1, MAP2K1, and LYN in the TCGA dataset. **p < 0.01; ***p < 0.001; ****p < 0.0001.

To validate the expression of two circRNAs (hsa_circ_0003763 and hsa_circ_0004928) and four immune-related mRNAs (FGF10, MET, RAC1, and LYN), blood samples from five GC patients and five normal controls were collected for qRT-PCR ( ). The clinical information of these individuals is listed in . The qRT-PCR result showed that RAC1 and LYN were significantly upregulated in GC, which was consistent with the integrated analysis results, while the expression trends of hsa_circ_0004928, MET, hsa_circ_0003763, and FGF10 were the same with the integrated analysis results without statistical significance.

Figure 13

In-vitro validation of hsa_circ_0003763, hsa_circ_0004928, FGF10, MET, RAC1, and LYN. **p-value < 0.01.

Table 6

Clinical information of individuals in the in-vitro experiment.

Group	Age (years)	Gender	Height (cm)	Weight (kg)	Pathological grading	Clinical staging			Differentiated degree	Metastasis	Metastatic site	Smoking history	Drinking history
					Stage I, II, III, IV	T	N	M	(G1, 2, 3, 4)
Control 1	63	Female	158	52	NA	NA	NA	NA	NA	NA	NA	No	No
Control 2	45	Female	156	54.5	NA	NA	NA	NA	NA	NA	NA	No	No
Control 3	66	Female	161	57	NA	NA	NA	NA	NA	NA	NA	No	No
Control 4	39	Female	159	59	NA	NA	NA	NA	NA	NA	NA	No	No
Control 5	59	Female	153	57	NA	NA	NA	NA	NA	NA	NA	No	No
GC 1	66	Male	170	73.2	IA	pT1a	0	0	G3	No	No	No	No
GC 2	58	Male	173	49	IIIB	pT4a	3a	0	G1	No	No	No	No
GC 3	70	Male	160	48	cIVB	cT4a	x	1	G3	Yes	Cervical lymph node	No	No
GC 4	58	Male	168	63.9	IV	cTx	x	1	Gx	Yes	Liver	Yes	Yes
GC 5	54	Female	158	56.6	IV	cT4	x	1	Gx	Yes	Abdominal cavity, left adnexa, and ascites	Yes	No

NA, not applicable.

Discussion

Up to now, no articles about hsa_circ_0050102 have been reported in any diseases. hsa_circ_0003763 was associated with pancreatic ductal adenocarcinoma (22). hsa_circ_0003763 was upregulated in hepatocellular carcinoma tissues and facilitated the invasion of hepatocellular carcinoma cells (23). hsa-miR-4537 was associated with vitreoretinal lymphoma (24). hsa-miR-145-5p played key roles in GC epithelial cells and affected cell invasion and migration of GC cells (25). In GC cells, hsa-miR-145-5p has been verified as the direct target of has_circ_0000376 (26). It has been demonstrated that hsa_circ_DLST can act as the sponge of has-miR-502-5p to regulate the NRAS/MAP kinase/ERK kinase 1/2 signaling pathways in GC cells (27). FGF10, a secreted factor, stimulated the proliferation of lung cancer cells (28). In GC, FGF10 was involved in several signaling pathways (29). Tgd cells were related to the prognosis of patients with stomach adenocarcinoma (30). In this study, we found regulatory pairs of hsa_circ_0050102–hsa-miR-4537–NRAS–Tgd cells and hsa_circ_0003763–hsa-miR-145-5p–FGF10–StromaScore in GC, among which, NRAS and FGF10 had potential diagnostic value for patients. In addition, NRAS and FGF10 were involved in both Rap1 and Ras signaling pathways in GC. Rap1 played important roles in the metastasis and invasion of various tumor cells by regulating cytoskeleton remodeling. In GC, some differentially expressed angiogenesis-related genes were associated with the Rap1 signaling pathway (31). Genes involved in the Ras signaling pathway were found in approximately 40% of GC patients (32). It has been suggested that regulatory networks of hsa_circ_0050102–hsa-miR-4537–NRAS–Tgd cells and hsa_circ_0003763–hsa-miR-145-5p–FGF10–StromaScore may play important roles in the proliferation, migration, and invasion of GC by involving signaling pathways of Rap1 and Ras.

hsa_circ_0040573 was found in human umbilical vein endothelial cell line (33). It was shown that RAC1 was an essential effector of GC malignant transformation and metastasis (34, 35). In gastric cancer tissues, the expression of RAC1 was increased, which was significantly related to TNM stage (36). Herein, we found the regulatory pairs of hsa_circ_0040573–hsa-miR-3686–RAC1–Tgd cells in GC. It was noted that RAC1 had potential diagnostic value for patients. Moreover, RAC1 was associated with Fc gamma R-mediated phagocytosis and cAMP signaling pathway. Fc gamma R-mediated phagocytosis was found in plasma exosomes from GC patients (37). It has been demonstrated that there was a significant positive correlation between nuclear factor kappa B (NF-κB) and cAMP-regulated phosphoprotein expression level in GC tissues (38). Thus, it can be seen that hsa_circ_0040573, RAC1, and Tgd cells played a crucial role in malignant transformation, metastasis, and TNM stage of GC.

The expression of hsa_circ_0001789 was found in GC (39). hsa-miR-1269b was downregulated in basal cell carcinoma and GC (40, 41). MET can be internalized by macrophages to educate them toward a protumorigenic phenotype (42). In addition, the upregulated paracrine hepatocyte growth factor can bind the c-MET receptor on the migrated GC cells to facilitate the proliferation of metastatic GC cells (43). In fat tissues, adipocytes regulated cancer development by their effects on the microenvironment. A previous report preliminarily showed that omental adipocytes could promote GC cell invasiveness (44). Herein, we found a regulatory pair of hsa_circ_0001789–hsa-miR-1269b–MET–adipocytes in GC. MET had a potential diagnostic value for patients. Furthermore, MET was involved in proteoglycans in cancer in GC. In gastric epithelia, the basal surface of cells was surrounded by the basement membrane, which was mainly composed of glycoproteins (45–47). The proteoglycan signaling pathway has been found in GC (48). It was indicated that the regulatory pair of hsa_circ_0001789–hsa-miR-1269b–MET–adipocytes played roles in cell proliferation and invasiveness of GC.

hsa_circ_0001013, upregulated in GC, regulated the expression of fibrillin 1 through competing with miRNA response elements of hsa-miRNA-182-5p, which led to metastasis in GC (49). hsa-miR-485-3p, a key regulator of gastrointestinal stromal tumors, served as a molecular biomarker and potential therapeutic target for this malignant disease (50). The role of hsa-miR-485-3p has been found in GC cells (51). MAP2K1 was involved in the PD-L1 pathway in GC (52). In our study, we found the regulatory network of hsa_circ_0001013–hsa-miR-485-3p–MAP2K1–Tgd cells in GC, among which, MAP2K1 had a potential diagnostic value for GC patients. In addition, MAP2K1 was enriched in T-cell receptor signaling pathway. Our result indicated that hsa_circ_0001013, hsa-miR-485-3p, MAP2K1, and Tgd cells may be associated with tumor metastasis in the development of GC.

hsa-miR-5584-3p played important roles in tongue squamous cell carcinoma by interacting with hsa_circ_087212 (53). LYN was significantly downregulated in gastric gastrointestinal stromal tumors with high-grade malignancy (54). In this study, we found the relationship between hsa_circ_0006089, hsa-miR-5584-3p, and LYN. Moreover, LYN had a potential diagnostic value and was negatively correlated with neurons. In addition, LYN was involved in chemokine signaling pathway. It was suggested that chemokines and their specific receptors can play an important role in GC progression via promotion of angiogenesis, invasion, survival, and metastasis (55). This indicated that the regulatory network of hsa_circ_0006089–hsa-miR-5584-3p–LYN–neurons played important roles in angiogenesis, invasion, survival, and metastasis of GC.

Besides the above circRNA–miRNA–immune-related mRNA regulatory network, some other circRNA–miRNA–mRNA regulatory networks were also found in GC, such as hsa_circ_0045602–hsa-miR-4538–WNT5A and hsa_circ_0018004–hsa-miR-199a-5p–KL. hsa_circ_0045602 was upregulated in GC tumor tissues (56). The expression of hsa-miR-4538 was decreased in acute myeloid leukemia (57). hsa-miR-199a-5p level was significantly increased in GC tissues. Moreover, a higher miR-199a-5p expression level of hsa-miR-199a-5p promoted the invasion of GC cells and was related to increased likelihood of lymph node metastasis (58). Thus, it can be seen that these molecules may be involved in tumor cell invasion and metastasis of GC.

In conclusion, our study found several ceRNA (circRNA–miRNA–immune-related mRNA) regulatory networks including hsa_circ_0050102–hsa-miR-4537–NRAS–Tgd cells, hsa_circ_0001013–hsa-miR-485-3p–MAP2K1–Tgd cells, hsa_circ_0003763–hsa-miR-145-5p–FGF10–StromaScore, hsa_circ_0001789–hsa-miR-1269b–MET–adipocytes, hsa_circ_0040573–hsa-miR-3686–RAC1–Tgd cells, and hsa_circ_0006089–hsa-miR-5584-3p–LYN–neurons in GC, among which FGF10, MET, NRAS, RAC1, MAP2K1, and LYN had a potential diagnostic value for GC patients. In addition, some signaling pathways were identified, such as Rap1 and Ras signaling pathways (involved NRAS and FGF10), Fc gamma R-mediated phagocytosis and cAMP signaling pathway (involved RAC1), proteoglycans in cancer (involved MET), T-cell receptor signaling pathway (involved MAP2K1), and chemokine signaling pathway (involved LYN). Our study may provide a novel field for understanding the molecular mechanisms of GC at the immunological levels. However, there are limitations to our study. Firstly, larger numbers of samples in qRT-PCR are needed. Secondly, an in-vitro gastric cell line experiment in cultures and tissues from healthy and gastric cancer patients is further needed to confirm the results of the in-silico analyses. Thirdly, the cell reporter assay for miRNA–circRNA–mRNA interaction is further needed.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

Ethics Statement

All participating individuals provided informed consent with the approval of the Ethics Committee of Zhejiang Hospital (2021-164K). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

GZ was the major contributor in the subject design. GZ contributed administrative support. PL provided the materials and samples and collected and sorted the data. ZW and GZ contributed to the data analysis and interpretation. All authors read and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.