Construct dysregulated miRNA-mRNA interaction networks to conjecture possible pathogenesis for Stomach adenocarcinomas

2.1Data collection from GEO and TCGA databases

To acquire candidate miRNA and mRNA datasets, we searched in the Gene Expression Omnibus (GEO) database with ‘(gastric OR stomach) AND (tumor OR cancer OR adenocarcinoma OR carcinoma) NOT (cell line)’ as the key word. Only microarray datasets measuring miRNA or mRNA expression levels in homo sapiens between normal tissue and tumor tissue were included. Under these constrains, totally 23 mRNA datasets and 8 miRNA datasets (GSE26595; GSE23739; GSE28700; GSE63121; GSE77380; GSE2564; GSE93415; GSE54397; GSE103236; GSE118897; GSE130823; GSE13861; GSE13911; GSE19826; GSE2685; GSE26899; GSE27342; GSE29272; GSE29998; GSE33335; GSE33651; GSE49051; GSE51575; GSE55696; GSE56807; GSE63089; GSE65801; GSE66229; GSE78523;GSE79973; GSE96668) listed in Table 1 were introduced for further analysis. RNA-sequencing profile and clinical information of gastric cancer collected in TCGA database were downloaded from the GDC data portal in the National Cancer Institute (https://portal. gdc.cancer.gov/). Besides, the overlapped differently expressed miRNAs involved in public databases including Human MicroRNA Disease Database (HMDD), Database of Differentially Expressed MiRNAs in human Cancers (dbDEMC) and miRCancer were also absorbed. The flowchart of the whole research is exhibited in Fig. 1.

2.2Data processing of RNA expression profile

GEO2R, an R-based online application using the GEOquery and limma R packages, was employed in differential analysis of GEO datasets. We performed Wilcoxon rank-sun test on sequencing data obtained from TCGA database with an R-based package called psych. Fold change (FC) and P-value were calculated in the two methods mentioned above. P-value < 0.05 and |log₂FC|> 0.58 were considered cut-off criterion for identifying differently expressed miRNAs (DE-miRNA) and mRNAs (DE-mRNA). Then, the intersection was taken to obtain the list of DE-miRNAs and DE-mRNAs in the GEO database and TCGA database. By summarizing GEO, TCGA, HMDD, dbDEMC, and miRCancer database, the list of DE-miRNAs and DE-mRNAs was ultimately identified.

2.3Functional enrichment analysis and speculation of agents

Gene oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DE-miRNAs were performed by DAVID-mirPath GO/KEGG analysis of DE-mRNAs were conducted on Hiplot (https://hiplot.com.cn) and presented in bubble plot P-value < 0.05 and counts > 2 were considered statistically significant. Connectivity MAP established connection between mRNA expression signatures and approximately 5,000 small molecule compounds, which helps to offer potential agents for gastric cancer treatment. Negative log10 transformed FDR q-values (fdr_q_nlog10) > 2 and normalized connectivity score (norm_cs) <-1 were set as the screening criterion.

2.4Predicting target mRNAs of DE-miRNAs

TarBase and MiRTarBase, databases collecting experimentally verified miRNA-mRNA interactions, were applied in assessing target mRNAs. Only miRNA-mRNA axes recorded in both databases were taken into consideration. To make sure that these negatively regulatory relationship between miRNA and mRNA can be observed in gastric tissue, we searched each miRNA-mRNA pair in PubMed. MiRNA-mRNA networks were plotted with Cytoscape v3.8.2.

2.5Sample collection, RNA isolation and reverse transcription quantitative PCR

We collected 30 paired Formalin-Fixed and Paraffin-Embedded (FFPE) samples of STAD tissue and adjacent normal stomach tissue. All patients had their stomach surgery in the First Affiliated Hospital of Nanjing Medical University None of the patients underwent preoperative treatment. Patients with stage IV gastric cancer do not have indications for surgery, so the samples collected in this study did not include stage IV patients. All patients signed written consent after fully understanding this research. This study was approved by the Institutional Review Boards of the First Affiliated Hospital of Nanjing Medical University (ID: 2016-SRFA-148). The clinical features of patients were disclosed in Table 2.

All stomach FFPE specimen were stored at 20^∘C. We isolated total RNA from FFPR stomach samples with RNAprep Pure FFPE Kit (TIANGEN Biotech, Beijing, China, DP439) according to the manual. Poly (A) Polymerase Kit (Takara) was applied in adding a poly (A) tail to each RNA. Reverse transcription was performed with PrimeScript RT reagent Kit (Takara, Tokyo, Japan, RR037A) at 37^∘C for 15 minutes, then 85^∘C for 5 seconds. Complementary DNA (cDNA) were stored at 4^∘C before conducting polymerase chain reaction (PCR). Quantitative PCR was carried out with SYBR Premix Ex Taq II (Takara, Tokyo, Japan, RR820A). We performed qPCR in a volume of 10 μL on the qTOWER384 (Analytik Jena) plate with the temperature set as 95^∘C for 20 s, followed by 40 cycles composed of 10 s at 95^∘C and 20 s at 60^∘C. The mRNA results were normalized with GAPDH and 18⁢S [17]. RNU6B was set as reference for miRNA [18]. We queried Primer Bank (https://pga.mgh.harvard.edu/primerbank/) for the primer sequence listed in Table 3.

Figure 2.

TarBase and MiRTarBase screened totally 18 negatively regulated miRNA-mRNA axes between 38 DE-miRNAs and 228 DE-mRNAs, which were plotted via Cytoscape. Orange nodes represent the upregulated miRNAs/mRNAs in STADs versus NCs, while blue nodes represent the downregulated miRNAs/mRNAs in STADs versus NCs.

Figure 3.

(A) The expression levels of KIT, LIFR, and ELL2 declined in STAD tissues. ANXA2 was overexpressed in tumor samples. (B) The miRNA expression level of miR-21-5p, miR-20a-5p, miR-221-3p, and miR-301a-3p were increased in STAD specimen. MiR-26a-5p, miR-29c-3p, and miR-1-3p were down regulated in STAD tissues. Data were presented as mean ± SEM. ^{*p <} 0.05, ^**p < 0.01, ^***p < 0.001 and ^****p < 0.0001 (Wilcoxon rank sum test). (C and D) The networks of miR-1-3p (down)/ANXA2 (up) and miR-301a-3p (up)/ELL2 (down) were the negatively regulated pairs verified via RT-qPCR.

Figure 4.

Receiver operating characteristic (ROC) curves of TCGA database and internal validation were separately presented in Fig. 4A and B. (A) The AUC of the complex model containing both miR-1-3p/ANXA2 and miR-301a-3p/ELL2 axes were 0.9208, the maximum among all 15 diagnosis models. (B) DCA presented great diagnostic performance of miRNA-mRNA axes in TCGA database and internal validation.

2.6Diagnostic efficiency assessment of MiRNA-mRNA axes and Survival analysis for STAD

We conducted binary logistic regression on miRNA and mRNA expression data with SPSS Statistics 26.0 (IBM, New York, USA). The receiver operation characteristic (ROC) curve and Decision Curve Analysis (DCA) were performed for comparing diagnostic efficiency between multiple models. ROC curves were plotted with data downloaded from TCGA database. According to survival data downloaded from TCGA database, we drew Kaplan-Meier survival curve with Hiplot. The average follow-up time was 577 days (1.58 years) with the maximal follow-up time equal to 3720 days (10.19 years). The medium divided all patients into low expression group and high expression group.

2.7Clinical subgroup analysis for STAD

Clinical information of stomach adenocarcinoma patients was downloaded from GDC data portal. According to their age (at first diagnosis), gender, tumor grade and tumor stage, cases were divided into several subgroups. Expression levels of miRNA-mRNA axes were compared between these subgroups by Wilcoxon rank-sum test. P-value < 0.05 was regarded as significant cut-off.

2.8Correlation analysis between miRNA-mRNA axes and tumor phenotypes

Wilcoxon Rank-sum test was applied in comparing the percentage of 22 immune cells between tumor samples and normal samples according to CIBERSORT (https://cibersort.stanford.edu/index.php/). Correlation analysis between proportion of discrepant immune cells and miRNA/mRNA expression level was exhibited in correlation heatmap plotted by Hiplot with Spearman’s correlation. ESTIMATE provided the Stromal Score, immune Score, estimate Score and tumor purity with mRNA expression signatures recorded in TCGA database. Tumor mutational burden (TMB) was measured by Masked Somatic Mutation data obtained with the ‘maftools’, an R-based package. We carried out Spearman’s correlation analysis between miRNA-mRNA axes and the tumor features mentioned above.

2.9Statistically analysis

MiRNA and mRNA expression data obtained from RT-qPCR was managed with 2 ^{- ΔΔCT} method (ΔCT = Ct_{miRNA or mRNA}- Ct_normalizer; Ct: the threshold cycle). The correlations between miRNA and mRNA expression data were calculated by Spearman’s correlation via psych package. GraphPad prism v8.0.1 was employed in ROC curves plotting and the area under the ROC curves (AUC) calculating. DCA curves were drawn with an R-based package named ‘rmda’. P-value < 0.05 was considered statically significant. Hiplot, Cytoscape, GraphPad prism and R language were applied in figure construction.

3.1Screening DE-miRNAs and DE-mRNAs in TCGA and GEO for STAD

Differential analysis of TCGA database provided basic DE-mRNAs and DE-miRNAs for STAD. As is presented in Table 1, a total of 23 mRNA GEO datasets and 8 miRNA GEO datasets met our inclusion criterion. Table 1 displayed the GEO accession, number of samples, researchers’ nationality, platform, and PMID of these datasets. 105 upregulated mRNAs and 123 downregulated mRNAs showed statistical significance in after calculating with GEO2R. Compared with miRNA expression profile of normal tissue, 22 upregulated and 16 downregulated miRNAs were considered significant. Combining with dbDEMC, HMDD and miRCancer, 22 upregulated and 16 downregulated miRNAs were selected ultimately.

3.2GO/KEGG analysis and potential compounds prediction

Outcome of function annotation and pathway enrichment analysis based on DE-mRNAs were presented in Fig. S1A and S1B. Mitotic cell cycle phase transition, extracellular matrix structural constituent and collagen-containing extracellular matrix were the GO terms of upregulated mRNAs with the most counts. Muscle system process and oxidoreductase activity acting on CH-OH group of donors were the most notable GO term of downregulated mRNAs. KEGG enrichment analysis suggested human papillomavirus infection crucial for upregulated mRNAs and gastric acid secretion remarkable for down regulated mRNAs. Outcomes of GO/KEGG analysis based on downregulated and upregulated miRNAs were separately revealed in Fig. S1C and S1D. Cellular amino acid metabolic process and transcription coactivator activity had smallest P-value for downregulated miRNAs enrichment analysis. KEGG analysis revealed that downregulated miRNAs were enriched in ErbB signaling pathway.

Following the DE-mRNAs signature query, small molecule agents including ixazomib, gemcitabine and tipifarnib were determined as potential therapeutic. The 15 compounds with the smallest connectivity scores were displayed in Fig. S2A. Agents’ mechanism of actions and their relative gene counts were exhibited in Fig. S2B.

Figure 5.

Patients were divided into different groups according to their clinical features, such as gender, age, tumor grade, and tumor stage. Subgroup analysis based on TCGA database suggested that ANXA2 expression was differed between patients older than 65 year and patients younger than 65. ELL2 was overexpressed in tissue with higher grade and stage.

Figure 6.

(A) This bar plot exhibited cell fraction of 22 immune cells in STAD versus normal tissue. (B) The relationship between differently infiltrated 8 immune cells and miRNA/mRNA axes was analyzed by Spearman’s correlation. (C) The association between expression levels of miRNA/mRNA axes and tumor microenvironment features including tumor mutation burden, tumor purity, ESTIMATE score, immune score and stromal score.

3.3Identification of miRNA-mRNA networks for STAD

To predict the role of DE-miRNAs in STAD, we collected their potential target mRNAs from TarBase and MiRTarBase. Figure 2 showed the miRNA-mRNA network. After searching in PubMed, totally 12 miRNA-mRNA axes were tested in gastric FFPE tissue via RT-qPCR.

3.4Corroboration of the DE-mRNAs and DE-miRNAs with RT-qPCR

Specific fold change and p-value of all miRNAs and mRNAs calculated via 2 ^{- ΔΔCt} method was listed in Table 4. MiRNA and mRNA expression data were exhibited in Fig. 3A and B separately. KIT (p< 0.0001), LIFR (p= 0.0007) and ELL2 (p= 0.0004) showed significant downregulated trend in gastric cancer tissue. Only ANXA2 (p= 0.0002) was upregulated in tumor tissue. Hsa-miR-26a-5p (p= 0.0005), hsa-miR-29c-3p (p= 0.029) and hsa-miR-1-3p (p= 0.0164) were observed downregulated in STAD compared with normal controls. Hsa-miR-21-5p (p< 0.0001), hsa-miR-20a-5p (p= 0.0175), hsa-miR-221-3p (p= 0.0467) and hsa-miR-301a-3p (p= 0.0324) expression level of gastric tumor tissue were remarkably higher than that of adjacent normal tissue. Spearman’s correlation analysis confirmed two reverse regulatory miRNA-mRNA axes. As is shown in Fig. 3C, hsa-miR-301a-3p was negatively regulated with ELL2 (p= 0.0123, r=-0.3212). Hsa-miR-1-3p expression was significantly correlated with ANXA2 (p< 0.0001, r=-0.5722) (Fig. 3D).

3.5Evaluating the diagnosis value of miRNA-mRNA axes in STAD

None of the GEO datasets contained both miRNA and mRNA expression profile of the same specimen. As a result, only expression patterns gained from TCGA data and RT-qPCR were applied in measuring predictive value. MiR-1-3p, ANXA2, MiR-301a-3p and ELL2 were combined as a panel with binary logistic regression ROC curves of 14 models and the AUC of each curve calculated with TCGA data were presented in Fig. 4A. The panel containing 4 signatures showed the best diagnosis value (AUC = 0.9208, 95% CI: 0.8804 to 0.9612). The DCA shown in Fig. 4B suggested satisfied diagnosis efficiency of the complex model in TCGA and RT-qPCR validation.

3.6Clinical subgroup analysis and survival analysis based on TCGA clinical information

GDC data portal provided clinical features and follow-up time of all clinical cases recorded in TCGA-STAD database. Patients were divided by age at first diagnosis, gender, tumor grade, and tumor stage. ANXA2 expression level was higher in patients with older first diagnosis age (p= 0.0269) (Fig. 5A). No significant difference of the miRNA-mRNA axes expression was observed between male and female patients (Fig. 5B). ELL2 expression was higher in patients with higher grade number (p= 0.0122) (Fig. 5C). ELL2 expression had a significant between early-stage (I + II) group and late-stage (III + IV) group (p= 0.0049) (Fig. 5D). Overall survival (OS) was employed to estimate tumor prognosis. As is demonstrated in Fig. S3, high ELL2 expression may suggest poor survival time (p= 0.041).

3.7Correlation analysis between miRNA-mRNA axes expression and tumor-related phenotypes

Among 22 types of immune cells, proportion of 8 immune cells showed significant difference between normal tissue and gastric tumor (Fig. 6A and Table 5). Results of Spearman’s correlation between immune cell proportion and miRNA/mRNA expression were presented in Fig. 6B. Percentages of M1 macrophages, monocytes, and activated memory CD4 T cells were correlated with miR-1-3p/ANXA2 axes. Regulatory T cells (Tregs) and follicular T helper cells were correlated with miR-301a-3p/ELL2 axes. Resting mast cells and M0 macrophages were significantly correlated with these 2 miRNA-mRNA axes mentioned above. ESTIMATE was applied to estimate the samples’ stroma and immunity levels. As is exhibited in Fig. 6C, the network of miR-1-3p and ANXA2 interacted with tumor mutational burden and tumor immune microenvironment. The p-value for Fig. 6B and C is presented in the supplementary table.