To identify GC data sets with relevant clinical information, systematic searches of GEO datasets (https://www.ncbi.nlm.nih.gov/geo/) were conducted. The inclusion criteria for the datasets were as follows: i) Sample size >50; ii) the studies presented relevant clinical information; and iii) datasets were generated using Affymetrix Human Genome 133 plus 2.0 Gene Chips (Affymetrix; Thermo Fisher Scientific, Inc.). In total, four data sets, GSE66229 (Cristescu et al, 2015) (n=400) (15), GSE15459 (Ooi et al, 2009) (n=200) (16), GSE57303 (Qian et al, 2014) (n=70) (17) and GSE34942 (Lei et al, 2012) (n=56) (18), were selected for further analysis. The batch function of SVA package (3.32.1 version) (https://bioconductor.org/packages/release/bioc/html/sva.html) was used to consolidate the 4 data sets. Data normalization was performed by R software (3.6.0 version; http://cran.r-project.org/) and gene expression levels were computed as mean values of all annotated probe sets (19).

The mRNA expression levels of FGF2 between tumor (n=626) and normal tissues (n=100) were compared. Both the median and the receiver operating characteristic (ROC) curve can be used as the grouping criteria. Due to the difference in the sample size between the groups, the median expression of FGF2 was selected to replace the ROC curve as the grouping standard. According to the median expression of FGF2, 612 patients, whose clinical information was available, were equally divided into two groups: An FGF2-high expression group (FGF2-H) and an FGF2-low expression group (FGF2-L). Differentially expressed genes [P<0.05, and |log fold change (FC)|≥1)] were identified between the two groups and analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using DAVID 6.8 (https://david.ncifcrf.gov/) (20). Clinical characteristics, including patient age, gender, tumor stage, Lauren classification (21), the extent of immune infiltration and overall survival (OS), were collected in order to identify correlations regarding FGF2 expression and clinical characteristics. The CIBERSORT (https://cibersort.stanford.edu/) method and the LM22 gene signature were used for immune infiltration analysis (22). All results were considered statistically significant at a threshold of P<0.05.

A Cox-prognosis model was constructed using R software (3.6.0 version). R packages including survival (version 2.44–1.1; http://github.com/therneau/survival), survminer (version 0.4.5; http://www.sthda.com/english/rpkgs/survminer/), survivalROC (contained in survcomp package; version 1.34.0; http://git.bioconductor.org/packages/survcomp) and survcomp (version 1.34.0; http://git.bioconductor.org/packages/survcomp) were used for the model and the calculation of C-index. In total, 612 patients with clinical information were randomly divided into two groups: A training group (n=306) and a validation group (n=306). The risk score for each patient was calculated based on the independent factors identified by the multivariate Cox model. The prognosis for each patient was evaluated according to risk scores. The sensitivity and specificity of the model were described using a time-dependent ROC curve (23).

Statistical analysis was performed using R software and SPSS version 23.0 (IBM Crop.). Graphical representations were generated using GraphPad Prism 7 (GraphPad Software, Inc.). χ² and Wilcoxon rank-sum tests were used for categorical and continuous variables, respectively. Multivariate Cox proportional-hazards analysis was used for identifying independent prognostic factors.

In total, four GEO datasets were extracted and normalized to analyze the expression patterns of FGF2 in GC. Compared with normal tissues, FGF2 mRNA expression levels were lower in GC (P<0.0001; Fig. 1). A total of 536 differentially expressed genes in the FGF2-H and FGF2-L cohorts were identified by R software and the results were analyzed using GO and KEGG (Fig. 2 and Table SI). GO analysis results revealed that the differentially expressed genes were predominantly involved in the regulation of ‘extracellular matrix organization’, ‘cell adhesion’ and ‘angiogenesis’ (Fig. 3A). KEGG analysis results revealed that the genes identified were predominantly involved in ‘focal adhesion’, ‘ECM-receptor interaction’ and the ‘cGMP-PKG signaling pathway’ (Fig. 3B). The biological functions and signaling pathways associated with these genes were closely related to the development of GC.

In total, 612 patients with clinical data were selected, and an association analysis of the FGF2 expression levels and clinical characteristics was performed. The patients were divided to the FGF2-H and FGF2-L groups, according to the median FGF2 expression. As listed in Table I, the patient age in the FGF2-L group was higher compared with that in the FGF2-H group. A statistical significance was found between stages II and IV (P<0.001; Table I), therefore, for further comparison, stages I and II were classified as early-stage and stages III and IV as advanced. The FGF2-H group had more patients with advanced-stage tumors compared with the FGF2-L group (P<0.001; Table I). Classifying the patient samples according to Lauren's criteria, a difference was observed in diffuse and intestinal types, but no difference was observed for the mixed types. In total, 366 patients were selected with detailed TNM staging in the four datasets for further analysis. There was a significant difference between the two FGF2-expressing groups in the T2 and T3 stages. Compared with FGF2-L, the number of patients with advanced-stage GC (T3 and T4) in the FGF2-H group was higher compared with the FGF2-L group (P<0.001; Table I). No significant difference was observed when comparing the N and M scores of the two FGF2-expressing groups (Table I).

Patients with P>0.05 were removed following CIBERSORT calculations (n=584; FGF2-H=301, FGF2-L=283), and non-parametric testing was used for comparison of the immune cell fractions (Figs. 4 and 5). The results revealed significant differences among the abundance of various immune cells, including CD4⁺ T cells, T regulatory cells, γδ T cells, macrophages, dendritic cells and mast cells (P<0.05).

Since FGF2 has a role in cell adhesion and because significant differences were identified in the abundance of different immune cell types in the tumor tissues, the survival of patients with GC may be affected by FGF2 expression levels. As shown in Fig. 6A, there was a strong association between higher FGF2 expression and shorter OS in patients with GC. Using multivariate analysis, FGF2 levels, age and tumor staging, were demonstrated to be independent prognostic indicators (Fig. 6B).

The risk score for each patient in the training group (n=306) and the validation group (n=306) was calculated based on the three independent factors, FGF2 levels, age and tumor staging, as identified by the multivariate Cox model (Fig. 6B). The survival analysis revealed that patients with low-risk scores had longer survival times in both the training and validation groups (P<0.001; Fig. 7A). ROC curves for the prognosis model were plotted, and the area under the curve was 78.2% in the training group and 82.8% in the validation group (Fig. 7B). The C-index was 0.725 in the training group [Hazard ratio (HR), 0.688–0.761; P<0.001] and 0.733 in the validation group (HR, 0.695–0.771; P<0.001). The results of the risk score and status plots (Fig. 7C and D) were consistent, and the accuracy and feasibility of the model were verified.