Analysis of microRNA expression profiles reveals a 5‑microRNA prognostic signature for predicting overall survival time in patients with gastric adenocarcinoma.

There is growing evidence supporting dysregulated microRNAs (miRNAs) as potential prognostic biomarkers in cancer. The present study aimed to identify an miRNA model set with prognostic power for patients with gastric adenocarcinoma. miRNA‑seq data from 155 patients and 37 controls were downloaded from The Cancer Genome Atlas (TCGA) database for a comprehensive analysis of miRNA expression profiles and were used as training data. A total of 5 prognostic miRNAs, which have not been previously reported, were identified using univariate and multivariate Cox regression analyses. A separate 155‑patient TCGA cohort was used as a validation set for evaluation of the risk model. Patients in the training set were assigned into high‑ and low‑risk groups according to the 5‑miRNA signature risk scores. Kaplan‑Meier survival analyses demonstrated that patients with high risk scores had significantly shorter survival times than those with low risk scores. The risk model validation confirmed the prognostic ability of this 5‑miRNA signature in predicting the risk status of patients. Stratification analysis for clinical prognostic variables demonstrated recurrence and age were significant prognostic factors in the low‑ and high‑risk groups, respectively. In conclusion, the present 5‑miRNA signature is a potential independent risk factor for patient outcomes. The risk model based on the 5‑miRNA signature performed well in predicting overall survival time in patients with gastric adenocarcinoma.

Introduction

Gastric cancer is one of the four most malignant tumors, accounting for ~10% of cancer-associated mortalities worldwide in 2015 (1–3). There are established treatments for gastric cancer, as well as therapies under development, yet it remains a lethal malignancy and was the second leading cause of cancer-associated mortalities in East Asia in 2012, due to the high morbidity and late diagnosis (2). The 5-year overall survival rate in gastric cancer, particularly the advanced and recurrent types, is <25% (4). Gastric or stomach adenocarcinoma (STAD) accounted for ~90% of gastric cancer cases worldwide in 2014 (5). The majority of patients with STAD in Western countries are diagnosed at advanced or metastatic stages (6). The early diagnosis of STAD or gastric cancer greatly improves the outcome of the patient. Therefore, diagnostic and prognostic biomarkers are urgently required for improving STAD diagnosis and predicting patient outcomes.

Prognostic markers implicate the close monitoring and treatment of high-risk patients to prolong their overall survival time (7). Traditional prognostic markers of gastric cancer used in clinical practice are principally clinicopathological variables, including age, tumor stage, Helicobacter pylori infection, response to chemotherapy and recurrence (8–11). The discoveries of novel prognostic biomarkers contribute to introducing and designing novel treatment strategies to improve patient survival. With the development and application of genetic engineering methods, large-scale genomic analyses have revealed various molecular signatures associated with gastric cancer outcomes, including gene mutations, mRNAs and non-coding RNAs [microRNAs (miRNAs/miRs) and long non-coding RNAs (lncRNAs)] (1,12–15). The association of a single molecule with a disease is limited to the complex mechanism of disease development, whereas multifactor signatures have exhibited superior diagnostic and prognostic abilities (7,13,15,16).

Over the past 5 years, various prognostic or diagnostic models based on sets of factors have been identified in gastric cancer, including several potential diagnostic miRNA signatures (7,13,15,16). Huang et al (1) identified a 6-miRNA signature comprised of 6 overexpressed miRNAs (miR-10b-5p, miR-132-3p, miR-185-5p, miR-195-5p, miR-20a-3p and miR-296-5p) detected in the serum of patients. Zhu et al (13) defined a 5-miRNA signature (miR-16, miR-25, miR-92a, miR-451 and miR-486-5p) as a potential diagnostic biomarker. The majority of these miRNAs, including miR-25, miR-92a, miR-132-3p, miR-296-5p, miR-195-5p, miR-451 and miR-486-5p, were associated with the development of gastric cancer and the survival time of the patients (17–23). The emergence of novel miRNA signatures with diagnostic and prognostic abilities has attracted a great amount of interest, and suggests that there is unlimited potential in mining valuable multifactor prognostic sets with predictive capacity for patients with solid tumors.

The present study was designed to identify a model prognostic miRNA set with predictive power in the outcome of patients with STAD. miRNAs associated with the prognosis of patients from The Cancer Genome Atlas (TCGA) database were identified using two-step Cox regression analysis. A predictive risk model based on the miRNA signature was established and validated using a sample-splitting method. The validation and assessment of the performance of the risk model was conducted using a Kaplan-Meier log-rank test and the area under the curve (AUC) following time-independent receiver operating characteristic (ROC) analysis. Stratification analyses were performed to assess the prognostic value of clinical variables. The potential of using the miRNA signature as a prognostic model for the outcome of patients with STAD was defined.

Materials and methods

Data collection

STAD miRNA-seq data, based on the Illumina HiSeq 2000 RNA Sequencing platform (Illumina, Inc., San Diego, CA, USA), were downloaded from TCGA database (https://portal.gdc.cancer.gov/) in May 10, 2018. Only data with information on patient survival and prognosis were selected (n=310), and the cases were randomly assigned into training and validation sets according to the analysis design (Fig. 1). Non-tumor samples (n=37) were assigned into the training group and were employed for the identification of differentially expressed miRNAs (DEmiRs), whereas the validation data were used for the evaluation and assessment of the risk model.

Figure 1.

Flow diagram of the data analysis process. TGCA, The Cancer Genome Atlas; STAD, stomach adenocarcinoma; miRNAs, microRNAs; DEmiRs, differentially expressed miRNAs; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification and hierarchical clustering analysis of DEmiRs

The DEmiRs between the STAD (n=155) and control samples (n=37) in the training set were identified using the edgeR package (version 3.20.9; http://bioconductor.org/packages/release/bioc/html/edgeR.html) (24). miRNAs were considered to be statistically significantly expressed in STAD samples with false discovery rate (FDR) <0.05 and |log2[fold-change(FC)]|≥0.5. DEmiRs were subjected to a two-way hierarchical clustering analysis using the centered Pearson's correlation algorithm in the Pheatmap package (version 1.0.8; http://cran.r-project.org/web/packages/pheatmap/index.html) (25,26). The analysis was performed in R (version 3.4.1; http://www.r-project.org/).

Selection of prognostic DEmiRs associated with the outcome of patients with STAD

Univariate and multivariate Cox regression analyses in the survival package (version 2.41.3; http://cran.r-project.org/web/packages/survival/index.html) (27) in R (version 3.4.1) were performed to define the prognostic DEmiRs. The hazard ratio (HR) and 95% confidence interval (CI) were estimated. The prognostic DEmiRs with a Kaplan-Meier log-rank test P-value of <0.05 were defined as independent prognostic factors for patients with STAD.

Establishment and evaluation of prognostic risk model

Step I: Determination of the optimal cut-off values of miRNA expression

Optimal cut-off values of the expression levels of the prognostic DEmiRs were defined using X-Tile Bio-Informatics software (version 2.41.3; http://medicine.yale.edu/lab/rimm/research/software.aspx) (28) based on the survival analysis (χ2 test). Monte Carlo sampling P<0.05 was set as the threshold for the cut-off value. The status of each DEmiR was defined as 0 (expression level < cut-off) or 1 (expression level > cut-off) (29).

Step II: Establishment of the risk model

The prognostic index, defined as the miR score or risk score of each sample, was calculated using the linear combination of the expression values weighted by the multivariate Cox regression coefficient (β) and expression status: miR score = Σβ miRNA n × status miRNA n, where status is 0 or 1 as previously defined, and n represents the miRNA name. The corresponding data were stratified into high- and low-risk groups according to whether their miR scores were higher or lower than the median.

Step III: Evaluation of the risk model

The prognostic difference between the high- and low-risk groups was analyzed using a Kaplan-Meier log-rank test, and the AUC of the time-independent ROC curve was used to evaluate the performance of the risk model in predicting high- and low-risk patients (30). P<0.05 was considered to indicate statistically a significant difference in the Kaplan-Meier log-rank test. The validation set was used for the evaluation and assessment for the performance of the risk model.

Analysis of risk factors

Univariate and multivariate Cox regression analyses were performed to define the independent prognostic risk factors for STAD, with the threshold of the Kaplan-Meier log-rank test being P<0.05. Stratification analyses of the potential clinical prognostic factors for patients with STAD were performed. Furthermore, Cox regression analyses were performed to identify the association between the clinical factors and the survival times of high- and low-risk patients. P<0.05 was considered to indicate statistically significant differences. A Kaplan-Meier survival analysis was performed for factors that were revealed to be significantly associated with the survival time of the patients.

miRNA-mRNA regulatory network and functional enrichment analysis

To identify the biological functions associated with the prognostic DEmiRs, a functional enrichment analysis was performed to identify the predicted targets of DEmiRs. Paired mRNA-seq data from patients assigned into the high- and low-risk groups were downloaded from TCGA, and differentially expressed genes (DEGs) were identified using the edgeR package with a cut-off of FDR<0.05 and |log2FC|≥0.5. Potential mRNA targets of the prognostic DEmiRs were predicted using TargetScan (version 7.2; http://www.targetscan.org/vert_72/) (31). Overlapping genes between the identified DEGs and the predicted targets of the DEmiRs were selected for the construction of an miRNA-mRNA regulatory network using Cytoscape (version 3.6.1; http://www.cytoscape.org/) (32). The DEG-associated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified using Gene Set Enrichment Analysis (GSEA; version 3.0; http://software.broadinstitute.org/gsea/index.jsp) (33), with P<0.05 considered to indicate statistical significance.

Results

Baseline characteristics of the patients

Data from a total of 310 patients with STAD and 37 healthy controls were included in the present study. Accordingly, data from 155 patients and the 37 controls were assigned as the training set and used for the identification of DEmiRs and definition of the risk model. The data from the remaining 155 patients were included into the validation set for validation of the risk model (Fig. 1). The baseline characteristics of the 310 patients are listed in Table I.

Table I.

Baseline characteristics of all patients with stomach adenocarcinoma.

Clinical characteristics	Training set (n=155)	Validation set (n=155)	Entire set (n=310)
Age, years (mean ± SD)	63.58±10.76	66.32±9.56	64.95±10.16
Sex, n
Male	99	103	202
Female	56	52	108
Reflux, n
Yes	20	16	36
No	77	85	162
N/A	58	54	112
Anti-reflux treatment, n
Yes	14	15	29
No	67	70	137
N/A	74	70	144
H. pylori infection, n
Yes	6	13	19
No	64	77	141
N/A	85	65	150
Radiation therapy, n
Yes	30	27	57
No	120	127	247
N/A	5	1	6
Metastasis stage, n
M0	142	142	284
M1	10	5	15
N/A	3	8	11
Node stage, n
N0	48	45	93
N1	41	41	82
N2	35	25	60
N3	29	40	69
N/A	2	4	6
Tumor stage
T1	8	18	26
T2	38	47	85
T3	61	78	139
T4	47	12	59
N/A	1	0	1
Pathological stage, n
I	25	7	32
II	49	26	75
III	67	76	143
IV	13	46	59
N/A	1	0	1
Grade, n
1	4	3	7
2	51	58	109
3	95	90	185
N/A	5	4	9
Recurrence, n
Yes	24	19	43
No	108	110	218
N/A	23	26	49
Survival, n
Succumbed	61	60	121
Alive	94	95	189
Overall survival time, months (mean ± SD)	18.74±17.36	18.56±17.53	18.65±17.45

SD, standard deviation.

Identification of DEmiRs

A total of 124 DEmiRs (Table SI) were identified following a comparative analysis of the miRNA-seq data from the training set (155 tumor samples and 37 controls) using the edgeR package, with the criteria of FDR<0.05 and |log2FC|≥0.5. The majority of the DEmiRs (88.71%; 110/124 miRNAs) were upregulated and 14 (11.29%) were downregulated (Fig. 2A). A two-way hierarchical clustering analysis of the DEmiRs revealed the distinct expression profiles of these miRNAs in tumor and control samples (Fig. 2B).

Figure 2.

Identification of DEmiRs and two-way hierarchical clustering analysis. (A) Scatter plot of the miRNAs from miRNA-seq data of patients with STAD. Red, green, and blue lines indicate the significantly upregulated, downregulated and non-significantly dysregulated miRNAs, respectively. DEmiRs were defined with a false discovery rate <0.05 and a |log2(fold-change)|≥0.5. (B) Heat map depicting the two-way hierarchical clustering of the 124 identified DEmiRs in the training set (155 STAD and 37 normal samples), analyzed using a centered Pearson's correlation algorithm. The color scale indicates the upregulation (red) or downregulation (green) of expression. STAD, stomach adenocarcinoma; UP, upregulated; DOWN, downregulated; NON, non-significantly regulated.

Identification of prognostic DEmiRs

The expression data of the 124 DEmiRs were subjected to a univariate Cox regression analysis with overall survival time as the dependent variable, and 13 potential prognostic DEmiRs were defined (log-rank test P<0.05; Table II). These potential prognostic DEmiRs were then subjected to a multivariate Cox regression analysis with overall survival time as the dependent variable, and 5 DEmiRs were ultimately identified to be prognostic miRNAs within the training group (P<0.05; Table II).

Table II.

Univariate and multivariate Cox regression analysis of miRNAs associated with survival of patients with stomach adenocarcinoma.

	Univariate analysis			Multivariate analysis
miRNA	β	HR (95% CI)	P-value	β	HR (95% CI)	P-value	Cut-off point[a]
hsa-mir-1255a	0.376	1.456 (0.941–2.252)	0.046[b]	0.827	2.286 (1.268–4.121)	0.006[b]	0.63
hsa-mir-3687	0.216	1.241 (0.966–1.593)	0.046[b]	0.360	1.433 (1.068–1.925)	0.017[b]	1.21
hsa-mir-548o	−0.520	0.595 (0.354–0.999)	0.025[b]	−0.800	0.449 (0.230–0.877)	0.019[b]	0.33
hsa-mir-9-3	0.390	1.477 (1.015–2.148)	0.021[b]	0.523	1.687 (1.083–2.629)	0.021[b]	0.68
hsa-mir-7-2	−0.274	0.761 (0.583–0.992)	0.022[b]	−0.381	0.683 (0.476–0.982)	0.039[b]	1.08
hsa-mir-216a	0.351	1.420 (1.061–1.901)	0.009[b]	0.342	1.408 (0.964–2.058)	0.077	N/A
hsa-mir-618	−0.520	1.339 (0.971–1.848)	0.038[b]	0.410	1.506 (0.938–2.418)	0.090	N/A
hsa-mir-504	0.306	1.357 (0.953–1.934)	0.046[b]	0.331	1.392 (0.882–2.197)	0.155	N/A
hsa-mir-556	−0.401	0.670 (0.455–0.983)	0.021[b]	−0.179	0.836 (0.522–1.339)	0.457	N/A
hsa-mir-493	0.268	1.308 (0.977–1.749)	0.036[b]	0.070	1.073 (0.684–1.683)	0.761	N/A
hsa-mir-1228	0.270	1.310 (0.951–1.803)	0.049[b]	0.042	1.043 (0.689–1.577)	0.844	N/A
hsa-mir-541	0.494	1.639 (0.996–2.695)	0.026[b]	−0.058	0.944 (0.502–1.774)	0.857	N/A
hsa-mir-496	0.409	1.506 (1.070–2.118)	0.009[b]	0.043	1.044 (0.609–1.787)	0.877	N/A

miRNA expression level cut-off point as determined by the X-tile Bio-Informatics software under the conditions of P<0.05.

P<0.05. miRNA, microRNA; HR, hazard ratio; CI, confidence interval; β, Cox regression coefficient.

The optimal cut-off values of the expression levels of the prognostic DEmiRs are listed in Table II and the graphical representation of the optimal cut-off points is displayed in Fig. 3. The expression of miRNAs hsa-mir-1255a (β, 0.827; HR, 2.286; 95% CI, 1.268–4.121; P=0.006), hsa-mir-3687 (β, 0.360; HR=1.433; 95% CI=1.068–1.925; P=0.017) and hsa-mir-9-3 (β, 0.523; HR, 1.687; 95% CI, 1.083–2.629; P=0.021) was negatively associated with the overall survival time; that of hsa-mir-548o (β, −0.800; HR, 0.449; 95% CI, 0.230–0.877; P=0.019) and hsa-mir-7-2 (β, −0.381; HR, 0.683; 95% CI, 0.476–0.982; P=0.039) was positively associated. Expressly, patients with high levels of hsa-mir-1255a, hsa-mir-3687 and hsa-mir-9-3 had high risk scores and short survival times, and patients with high levels of hsa-mir-548o and hsa-mir-7-2 had low risk scores and longer survival times (Fig. 3). These findings indicate that elevated expression levels of hsa-mir-1255a, hsa-mir-3687 and hsa-mir-9-3, combined with low levels of hsa-mir-548o and hsa-mir-7-2, are associated with a poor prognosis.

Figure 3.

Determination of the optimal cut-off values of 5 miRNAs using X-Tile Bio-Informatics software. Red and green X-tile plots indicate that high and low miRNA expression level, respectively, is associated with longer survival time. Larger low and high population indicates larger population with low and high miRNA expression level, respectively. Red and gray lines in the Kaplan-Meier survival plots represent the analysis of patients with low and high miRNA expression levels, respectively. The optimal cut-off information is listed in the box below the survival plot. The data displayed are of (A) miR-7-2, (B) miR-9-3, (C) miR-548o. (D) miR-1255a and (E) miR-3687. miRNA, microRNA; No, patient number.

Risk model training

The risk score was calculated with coefficients from the multivariate Cox regression analysis by incorporating the 5 prognostic miRNAs. The predictive risk model for the training data using the 5-miRNA signature was constructed using the formula: miR=(0.827) × Status_hsa-mir-1255a + (0.360) × Status_hsa-mir-3687 + (−0.800) × Status_hsa-mir-548o + (0.523) × Status_hsa-mir-9-3 + (−0.381) × Status_hsa-mir-7-2. Overall, 77 (49.68%) and 78 samples (50.32%) from the training set were assigned into low-risk (−1.18047<miR≤-0.354) and high-risk (−0.354<miR≤1.710) groups, respectively, according to the distribution and the median value (−0.354) of the risk score for the patients (Fig. 4A). The overall survival rate within the low-risk group during the follow-up period was 79.22% (61/77 patients), whereas that within the high-risk group was 42.31% (33/78 patients). The Kaplan-Meier log-rank test revealed that the patients in the high-risk group were associated with significantly shorter survival times than those in the low-risk group (HR, 3.347; 95% CI, 1.903–5.886; P<0.01; Fig. 4B).

Figure 4.

(A) Risk score distribution and (B) Kaplan-Meier survival analysis of the training set of patients with stomach adenocarcinoma. (C) The risk score distribution and (D) Kaplan-Meier survival analysis of the validation set of patients. The risk model defined the patients into high- and low-risk groups. A log-rank test was used to analyze the difference between the survival times of the patients from the high- (red) and low-risk (black) groups. Significant differences between the two groups were observed. HR, hazard ratio; miR score, risk score.

Risk model validation

Data from TCGA validation set were used to evaluate the performance and prognostic power of the predictive risk model based on the 5-miRNA signature. Fig. 4C displays the risk score distribution for patients with STAD based on the 5 prognostic miRNAs. Overall, 77 (49.7%) and 78 (50.3%) patients were assigned into the low-risk (−1.180<miR≤-0.021) and high-risk (−0.021<miR≤1.350) groups, respectively. During the follow-up period, the overall survival rate within the low-risk group was 70.1% (54/77 patients), whereas that within the high-risk group was 52.6% (41/78 patients). As expected, the Kaplan-Meier log-rank test demonstrated that the high-risk group exhibited shorter survival times. A significant difference in survival times was observed between the low- and high-risk groups (HR, 2.360; 95% CI, 1.39–4.008; P<0.01; Fig. 4D).

Assessment of the risk model

The performance of the 5-miRNA signature risk model was evaluated by constructing a time-independent ROC curve. Fig. 5 displays the ROC curves and AUCs of the risk model. The AUC for the training and validation sets was 0.939 and 0.901, respectively (Fig. 5A and B). These results indicate the high performance and prognostic ability of the 5-miRNA signature risk model in predicting high- and low-risk patients with STAD.

Figure 5.

Time-dependent receiver operating characteristic curves of the sensitivity and specificity of the 5-microRNA signature as a predictive marker for the risk status of patients with stomach adenocarcinoma from (A) the training set and (B) the validation set. AUC, area under the curve.

Validation of the entire cohort

A similar risk stratification was revealed when the 5-miRNA signature was applied to the entire TCGA STAD cohort (Fig. 6). The training and validation cohorts were pooled together and all patients were assigned into low- (n=155) and high-risk (n=155) groups, based on the 5-miRNA signature risk score of the patients. The Kaplan-Meier log-rank test revealed that the patients in the low-risk group had significantly longer survival times than those in the high-risk group (HR, 2.840; 95% CI, 1.937–4.162; P<0.01; Fig. 6A). The AUC of risk model for the entire cohort was 0.918 (Fig. 6B), demonstrating the predictive power of the 5-miRNA signature prognostic model in estimating the overall survival time of patients with STAD.

Figure 6.

Risk model of entire TCGA cohort with stomach adenocarcinoma. (A) Kaplan-Meier survival analysis of patients in the high-risk (red) and low-risk (black) groups. A significant difference was observed between the groups (P<0.01). (B) A time-independent receiver operating characteristic curve demonstrating the performance of the 5-microRNA signature as a predictive marker for the risk status of patients with stomach adenocarcinoma. HR, hazard ratio; AUC, area under the curve.

Prognostic value of clinical variables

The prognostic values of clinical variables were analyzed using univariate and multivariate Cox regression analyses (Table III). Overall, 7 of these parameters, including age, radiotherapy, pathological differentiation, classification, stage and recurrence, were identified in the univariate analysis. The multivariate analysis ultimately identified 3 independent prognostic clinical variables as risk factors, including age (HR, 1.953; 95% CI, 1.176–3.242; P=0.010), recurrence (HR, 2.529; 95% CI, 1.491–4.289; P=0.006) and radiotherapy (HR, 0.410; 95% CI, 0.198–0.849; P=0.016) (Table III; Fig. 7).

Table III.

Analysis of the prognostic value of clinical variables in the entire tested dataset (n=310) of patients with stomach adenocarcinoma.

	Univariate analysis		Multivariate analysis
Variables	HR (95% CI)	P-value	HR (95% CI)	P-value
Age, years (≤65/>65)	1.551 (1.077–2.234)	0.018[a]	1.953 (1.176–3.242)	0.010[a]
Sex (male/female)	1.462 (0.979–2.182)	0.061	N/A	N/A
Reflux (yes/no)	0.728 (0.375–1.412)	0.345	N/A	N/A
Anti-reflux treatment (yes/no)	0.842 (0.459–1.543)	0.577	N/A	N/A
H. pylori infection (yes/no)	0.435 (0.173–1.092)	0.069	N/A	N/A
Radiation therapy (yes/no)	0.469 (0.281–0.786)	0.003[a]	0.410 (0.198–0.849)	0.016[a]
Metastasis stage (M0/M1)	2.497 (1.301–4.795)	0.004[a]	2.569 (1.115–5.915)	0.267
Node stage (N0-N1/N2-N3)	1.605 (1.116–2.308)	0.010[a]	0.995 (0.507–1.953)	0.989
Tumor stage (T1-T2/T3-T4)	1.661 (1.050–2.627)	0.028[a]	1.703 (0.889–3.259)	0.108
Pathological stage (I–II/III–IV)	1.959 (1.329–2.887)	0.001[a]	1.548 (0.723–3.316)	0.261
Grade (1/2/3)	1.323 (0.906–1.933)	0.146	N/A	N/A
Recurrence (yes/no/N/A)	2.198 (1.358–3.557)	0.001[a]	2.529 (1.491–4.289)	0.001[a]

P<0.05. HR, hazard ratio; CI, confidence interval.

Figure 7.

Three independent clinical prognostic variables in patients with STAD. Kaplan-Meier survival curves according to (A) age, an independent risk factor, (B) radiotherapy, an independent protective factor, and (C) recurrence, an independent risk factor, for patients with STAD. Significant differences were observed between the groups for all three factors (P<0.05). STAD, stomach adenocarcinoma; HR, hazard ratio.

Stratification analysis of the prognostic clinical variables

Stratification analyses were performed for the 3 independent clinical prognostic variables in the high- and low-risk groups. The information regarding tumor-node-metastasis stage, neoplasm histological grade and pathological stage of the tumors was obtained from TCGA. All patients were assigned into high- (n=155) and low-risk (n=155) groups, according to their 5-miRNA signature risk score. In the univariate analysis, radiotherapy, pathological stage and recurrence were significant risk factors in the low-risk group; whereas age, pathological node status and pathological stage were significant in the high-risk group (P<0.05; Table IV). In the multivariate analysis, recurrence (HR, 3.852; 95% CI, 1.439–10.313; P=0.007) and age (HR, 1.696; 95% CI, 1.081–2.660; P=0.022) were the only significant prognostic variables in the low- and high-risk groups, respectively.

Table IV.

Stratification analysis of the prognostic value of clinical variables in low-risk and high-risk group.

	Low risk group (n=155)				High risk group (n=155)
	Univariate analysis		Multivariate analysis		Univariate analysis		Multivariate analysis
Variables	HR (95% CI)	P-value	HR (95% CI)	P-value	HR (95% CI)	P-value	HR (95% CI)	P-value
Age, years (≤65/>65)	1.936 (0.980–3.824)	0.053	N/A	N/A	1.609 (1.037–2.496)	0.032[a]	1.696 (1.081–2.660)	0.022[a]
Sex (male/female)	1.704 (0.806–3.599)	0.158	N/A	N/A	1.331 (0.826–2.144)	0.239	N/A	N/A
Reflux (yes/no)	0.806 (0.241–2.692)	0.725	N/A	N/A	0.703 (0.314–1.574)	0.389	N/A	N/A
Anti-reflux treatment (yes/no)	1.324 (0.503–3.484)	0.568	N/A	N/A	0.475 (0.198–1.141)	0.089	N/A	N/A
H. pylori infection (yes/no)	0.487 (0.145–1.638)	0.236	N/A	N/A	0.476 (0.114–1.986)	0.297	N/A	N/A
Radiation therapy (yes/no)	0.388 (0.165–0.912)	0.025[a]	0.455 (0.123–1.676)	0.236	0.639 (0.336–1.217)	0.169	N/A	N/A
Metastasis stage (M0/M1)	4.707 (1.401–5.821)	0.006[a]	3.808 (0.785–18.466)	0.397	1.847 (0.842–4.052)	0.126	N/A	N/A
Node stage (N0-N1/N2-N3)	1.570 (0.831–2.965)	0.161	N/A	N/A	1.676 (1.071–2.621)	0.022[a]	1.200 (0.675–2.135)	0.534
Tumor stage (T1-T2/T3-T4)	1.670 (0.736–3.789)	0.215	N/A	N/A	1.736 (0.965–3.122)	0.062	N/A	N/A
Pathological stage (I–II/III–IV)	2.136(1.063–4.293)	0.029[a]	2.106 (0.686–6.460)	0.193	1.960 (1.219–3.150)	0.005[a]	1.699 (0.917–3.147)	0.092
Grade (1/2/3)	1.087(0.577–2.047)	0.797	N/A	N/A	1.374 (0.848–2.226)	0.195	N/A	N/A
Recurrence (yes/no/N/A)	4.799 (1.924–11.97)	0.000[a]	3.852 (1.439–10.313)	0.007[a]	1.397 (0.784–2.49)	0.255	N/A	N/A

P<0.05. HR, hazard ratio; CI, confidence interval.

The Kaplan-Meier survival analyses revealed that patients without recurrence in the low-risk group survived longer than those with recurrence (HR, 4.799; 95% CI, 1.924–11.970; P<0.01; Fig. 8A), and patients aged >65 years in the high-risk group had notably shorter survival times than those ≤65 years old (HR, 1.609; 95% CI, 1.037–2.196; P=0.032; Fig. 8B). Together, the findings suggest that age and recurrence are independent risk factors for patients with STAD with higher and lower risk scores, respectively.

Figure 8.

Stratification analysis of the prognostic value of (A) recurrence and (B) age in the low- and high-risk patients, respectively. Kaplan-Meier survival analysis revealed significant differences between the groups with and without recurrence in the low-risk patients, and between younger and older patients in the high-risk group (P<0.05). HR, hazard ratio.

Functional characteristics of the prognostic miRNAs

Since the 5-miRNA signature was calculated to be an independent risk factor for patients with STAD, a functional analysis was performed on the targets of the 5 prognostic DEmiRs. The corresponding mRNA-seq data from patients in the high- and low-risk groups were downloaded and a total of 244 DEGs (FDR<0.05 and |log2FC|≥0.5; Fig. 9A) were identified, as were 86 predicted targets of the 5 prognostic DEmiRs using TargetScan. The 244 DEGs are listed in Table SII. The miRNA-mRNA regulatory network was comprised of 91 nodes and 119 interactions (Fig. 9B). The GSEA KEGG pathway enrichment analysis revealed 4 pathways involving target genes, including ‘T cell receptor signaling pathway’ (P=0.011), ‘Cytokine cytokine receptor interaction’ (P=0.011), ‘Cell adhesion molecules (CAMs)’ (P=0.021) and ‘Toll-like receptor signaling pathway’ (P=0.029) (Table V), indicating the potential roles of the 5 prognostic DEmiRs and their targets in the development of STAD.

Figure 9.

Differentially expressed genes between high- and low-risk groups, and the miRNA-mRNA regulatory network. (A) The DEGs from mRNA-seq data of patients in high-risk and low-risk groups were identified (FDR<0.05; |log2FC|≥0.5). Blue and red triangles indicate significantly downregulated and upregulated genes, respectively. (B) miRNAs and DEGs are indicated with squares and cycles, respectively. The color represents the log2FC values, as denoted by the scale. The mRNA targets of the 5 prognostic miRNAs were predicted using TargetScan version 7.2. Cytoscape version 3.6.1 was used to visualize the network. DEGs, differentially expressed genes; miRNA, microRNA; FDR, false discovery rate; FC, fold-change.

Table V.

Significant GSEA KEGG pathways associated with target genes of the 5 prognostic microRNAs.

Pathway	ES	NES	NOM P-value	Genes
‘T cell receptor signaling pathway’	0.508	1.692	0.011	ITK, PIK3CD, CD8B, CD8A, CD3G, CD3E, CD3D, RASGRP1,ZAP70, CARD11
‘Cytokine cytokine receptor interaction’	−0.399	−1.710	0.011	CXCL6, CSF3, CXCL5, BMP2, IL6, IL1B
‘Cell adhesion molecules’	0.561	1.537	0.021	ITGAL, CD2, CD6, ITGB7, CD8B, CD8A
‘Toll like receptor signaling pathway’	−0.646	−1.503	0.029	IL6, IL1B

KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis; ES, enrichment score; NES, normalized enrichment score; NOM, nominal.

Discussion

Outcome-associated molecular signatures have implications on prognosis and the molecular mechanisms underlying diseases and cancer types (34–37). In the present study, a comprehensive analysis of miRNA expression profiles was performed on TCGA datasets from patients with STAD. A 5-miRNA prognostic signature with prediction power for survival was identified using Cox regression analysis and sample splitting techniques. Of the 5 miRNAs, 3 (hsa-mir-1255a, hsa-mir-3687 and hsa-mir-9-3) were associated with short survival times, and 2 (hsa-mir-548o and hsa-mir-7-2) were associated with longer survival times in patients with STAD. The performance and prognostic ability of the 5-miRNA risk model was determined using a 155-sample TCGA validation cohort. Furthermore, the 5-miRNA signature, patient age and tumor recurrence were revealed to be independent risk factors for patients with STAD.

It has been reported that the incidence of gastric cancer rises progressively with age. The risk and occurrence of gastric cancer are low in individuals <30 years old, and they gradually increase with age, peaking at >50 years old (5,10,11,38,39). Studies have reported that older age (>60 years) and H. pylori infection are independent risk factors for gastric cancer, and that H. pylori infection is more prevalent in individuals of older ages (10,11). Several studies have reported that the eradication of H. pylori infection led to a lower incidence and recurrence of gastric cancer (40–42). It has been demonstrated that recurrence of malignancies affects the prognosis of patients (11,43). Local recurrence of malignant tumors may be associated with distant tumor metastasis, poor overall survival time and mortality (44–46). In the present study, the stratification analysis of patients with high and low risk scores suggested that age and recurrence were independent risk factors of a poor overall survival time in patients with high and low 5-miRNA signature risk scores, respectively. These findings support the high performance and prognostic power of the 5-miRNA signature risk model in predicting the overall survival time of patients with STAD.

Numerous oncogenic and tumor-suppressor genes, miRNAs, lncRNAs and individual signatures have been identified and demonstrated to be diagnostic or prognostic markers for patients with many types of cancer, including gastric cancer (13,16,47,48). In the present study, a 5-miRNA signature (hsa-mir-1255a, hsa-mir-3687, hsa-mir-9-3, hsa-mir-548o and hsa-mir-7-2) was identified as an independent prognostic predictor of survival time in patients with STAD. Among these 5 miRNAs, downregulated miR-9-3 in human colorectal cancer (CRC) (49), hsa-mir-548o in glioblastoma (50) and hsa-mir-7-2 in thyroid cancer (51) had been identified, as well as upregulated hsa-mir-1255a in cirrhotic hepatocellular carcinoma (52), hsa-mir-3687 in prostate cancer cell lines (53) and hsa-mir-7-2 in renal cell carcinoma (54). To the best of our knowledge, no study has reported either the association of these miRNAs with the prognosis of patients with cancer, or their association with STAD. The performance of the 5-miRNA risk model in predicting the survival time of high-risk patients suggests that the 5-miRNA signature is a novel prognostic marker in STAD.

In the prediction of the targets of the identified 5 miRNAs, Kelch repeat and BTB domain-containing protein 11 (KBTBD11) and calcium/calmodulin-dependent protein kinase type 1D (CAMK1D) were two common targets of hsa-mir-1255a and hsa-mir-3687 (Fig. 9). KBTBD11 has been reported to function as a putative tumor suppressor in CRC, with its knockdown leading to enhanced CRC cell proliferation (55). Gong et al (55) demonstrated that the KBTBD11 polymorphism rs11777210 regulates the binding with Myc proto-oncogene protein, a transcription factor that negatively regulates the expression of KBTBD11 (56). One study reported the elevation of CAMK1D in metastatic breast cancer, with its upregulation in breast epithelial cells triggering cell proliferation, epithelial-mesenchymal transition, migration and invasion abilities, and reduced cell adhesion (57). The upregulation of CAMK1D in gastric adenocarcinoma has been reported following gene expression analysis using a DNA microarray (58). These studies imply the oncogenic potential of KBTBD11 inhibition and CAMK1D expression. In addition, Fussek et al (53) revealed that miR-3687 was involved in cell cycle regulation and was elevated in the G0/G1 phase. In the present study, the upregulation of the KBTBD11 and CAMK1D genes was demonstrated in STAD samples (fold-change >0.7). Furthermore, CAMK1D was confirmed to be downregulated by the protective miRNA hsa-mir-548o (Fig. 9). These results indicate the complex and important roles of hsa-mir-548o, hsa-mir-1255a and hsa-mir-3687 in STAD development.

Krüppel-like factor 12 (KLF12) is a transcription repressor and a known participating factor in the progression of human gastric cancer (59). Nakamura et al (59) reported that KLF12 mRNA levels were associated with tumor size and progression of gastric cancer. The study demonstrated that KLF12 enhanced gastric cancer cell proliferation and invasion, and the selective knockdown of KLF12 in gastric cancer HGC27 cells, resulted into marked proliferation arrest by deregulating the expression of proliferation-associated genes. miR-137 is frequently inhibited in gastric cancer (60). It has been reported that miR-137 expression inhibits cell proliferation and migration and arrests cell cycle at the G0/G1 phase in gastric cancer cells by regulating KLF12 (61). Furthermore, Mak et al (62) demonstrated that KLF12 was regulated by miR-141, and the knockdown of KLF12 promoted cell proliferation, tumor growth, metastasis and anoikis resistance in ovarian cancer cells. In the present study, KLF12 was revealed to be regulated by two miRNA risk factors (hsa-mir-1255a and hsa-mir-9-3) and two protective miRNAs (hsa-mir-548o and hsa-mir-7-2). The target mRNAs of these 5 miRNAs were associated with the ‘Cell adhesion molecules (CAMs)’ pathway. These findings imply the complex mechanisms underlying KLF12-mediated cell proliferation, which may contribute to STAD prognosis.

In conclusion, the present study identified a novel 5-miRNA signature risk model (hsa-mir-1255a, hsa-mir-3687, hsa-mir-9-3, hsa-mir-548o and hsa-mir-7-2) with prominent performance in predicting high risk scores and the overall survival time of patients with STAD. These circulating miRNAs may be monitored as risk factors (hsa-mir-1255a, hsa-mir-3687, and hsa-mir-9-3) or protective indicators (hsa-mir-548o and hsa-mir-7-2) with prognostic ability in patients with STAD. Future efforts should focus on uncovering the molecular mechanism associated with the 5-miRNA signature in STAD.