The GEO database (http://www.ncbi.nlm.nih.gov/geo) was used to retrieve data with “stomach adenocarcinoma” as the keywords and human as the species. Datasets that covered cancer tissue and normal adjacent tissue, came from the same platform, and contained at least 20 samples were selected, and then, the gene expression profiles and their clinical data were downloaded. From TCGA portal (https://tcga-data.nci.nih.gov/tcga/), we collected the STAD RNA-seq data and related clinical parameters.

Gene level expression data were normalized and then log2 transformation is provided by the limma package of R software (version 3.6.3). For GEO datasets, data were analyzed using the GEO2R analysis tool, and the DEGs were identified at adjusted p < 0.05 and |Log₂FC| >1. For TCGA-STAD cohort, DEGs were identified with false discovery rate (FDR) < 0.05 and |Log₂FC| >1 via the edge R package (9).

Overlapping DEGs were screened based on the p-value and fold change (FC)/log(FC), and the top 50 genes were selected for Kaplan–Meier survival analysis. Log-rank p-values for Kaplan–Meier plots were calculated using an R package for survival analysis. Then, we screened genes with log-rank p < 0.05 as prognostic-related genes for subsequent analysis.

To establish the prognostic model of STAD, univariable and multivariable Cox regression analyses were performed on the prognostic-related genes by the survival R package (10). Owing to the lack of survival information on GEO, we randomly divided the patients with complete survival information in TCGA-STAD dataset into training set and validation set, fit the model in the training set, and assessed its performance in the validation set. Then, the prognostic model was established based on corresponding coefficients of the prognostic genes of STAD.

Further, the training set was divided into high- and low-risk groups according to the median value of risk score. Kaplan–Meier survival analysis was performed to estimate overall survival (OS) between the two groups by the survival R package. Time-dependent receiver operating characteristic (ROC) curves were plotted by time ROC R package, and the areas under ROC curves (AUCs) were calculated to test the efficiency of the prognostic model (11). Univariable and multivariable Cox regression analyses were applied on clinical data (including age, gender, TMN stage, and histologic grade) and risk scores to assess whether the risk model was an independent prognostic factor of clinical parameters.

Interaction network analysis was obtained by employing STRING v11.5 database (http://string-db.org/), keeping default parameters. The topological properties of the PPI network included average shortest path length, betweenness centrality, closeness centrality, degree, eccentricity, neighborhood connectivity, radiality, stress, and topological coefficient. Molecular complex detection (MCODE) analysis was applied to the prognostic-related gene network to identify densely connected subnetwork modules. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed to identify significant pathways via the “cluster Profiler” package in R (12). The items of biological processes were further analyzed by GO classifications. The adjusted p < 0.05 was considered to indicate a statistically significant difference. In addition, gene set enrichment analysis (GSEA) was utilized to determine the core gene-related signaling pathways by the “cluster Profiler” package in R. Results with absolute value of normalized enrichment score > 1, FDR < 0.25, and adjust p < 0.05 were considered statistically significant. 1D linear domain structures and 3D structures of proteins were visualized using cBioPortal (http://www.cbioportal.org/).

Differential expression of the prognostic genes between normal and STAD-related tissues was verified. Additionally, immunohistochemistry staining of the prognostic genes in STAD and normal tissues was acquired from the Human Protein Atlas database (https://www.proteinatlas.org/). According to the results of univariate and multivariate Cox regression analyses, a nomogram was created using the rms and survival package of R (13). Additionally, a calibration plots were generated to validate the nomogram.

The clinical data of GSE27342, GSE63089, and TCGA-STAD were shown in Table 1. We found 474 DEGs in GSE27342 profile (287 upregulated and 187 downregulated, Figure 2A), 732 DEGs in GSE63089 profile (622 upregulated and 110 downregulated, Figure 2B), and 5,494 DEGs in TCGA-STAD cohort (2,659 upregulated and 2,835 downregulated, Figure 2C). Subsequently, a total of 229 overlapping DEGs (159 upregulated and 70 downregulated) were screened among the three datasets (Figure 2D and Additional File 1).

We selected the top 50 overlapping DEGs (cutoff: p < 0.05 and |Log₂FC| >1.75) as candidate genes. According to log-rank p < 0.05 by Kaplan–Meier survival analysis, 11 genes (ADAM2, BGN, COL10A1, MMP1, MMP7, MFAP2, CTHRC1, P4HA3, SFRP4, TNFRSF11B, and FAP) were screened as prognostic-related genes for following research, and the Kaplan–Meier plots were shown in Figure 3.

Since the expression value of ADAM2 was zero in half of the samples, it was impossible to group by the median. The results of the univariate and multivariate proportional hazards regression analyses of the 10 prognostic-related genes associated with clinical outcomes are shown in Table 2. Multivariate regression analysis revealed COL10A1, MFAP2, CTHRC1, P4HA3, and FAP as the risk factor, and the risk score formula for OS was as follows:

In addition, the patients of TCGA-STAD dataset were divided into training set and validation set. Training set consisted of 186 STAD cases whereas validation set consisted of 185 STAD cases. Patients with STAD were divided into high- and low-risk subgroups according to the median value of risk score (cutoff = 14.9). In the training set, the survival analysis showed that the OS rates in the high-risk group were significantly lower than those in the low-risk group (p = 0.0025, Figure 4C). The time-dependent ROC curves offered a survival prediction that the AUCs were 0.576 (1-year OS), 0.733 (3-year OS), and 0.887 (5-year OS). Result showed that the risk model had a good ability to predict long-term prognosis of STAD (Figure 4B). The heatmap showed that the expression levels of five core genes were higher in patients with STAD with high-risk scores than those with low-risk scores (Figures 4A, 5A). Meanwhile, data in the validation set showed the similar results: OS rates in the high-risk group were significantly lower than those in the low-risk group (p = 0.045, Figure 5C); the time-dependent ROC curves (Figure 5B) predicted that the AUCs were 0.530 (1-year OS), 0.599 (3-year OS), and 0.702 (5-year OS). Moreover, using multivariate Cox regression analysis, the prognostic model was identified as an independent predictor for patients with STAD (p = 0.008, Table 3).

The PPI network of the five core genes was shown in Figure 6A. The topological properties of the PPI network for each gene were shown in Additional File 6. Highly interconnected subcluster of the five core genes was shown in Figure 6B, and the subcluster consisted of 53 nodes and 305 edges (score = 11.731), which represented relatively stable of protein in the network.

The results of GO enrichment analysis (Figure 7A) showed that the five core genes significantly focused on extracellular matrix organization, extracellular structure organization (biological process, BP); extracellular matrix structural constituent, dipeptidyl-peptidase activity (molecular function, MF); and collagen-containing extracellular matrix, collagen trimer (cell components, CC). Meanwhile, we found that in terms of biological processes, the genes were mainly focused on extracellular matrix, cell cycle, and Wnt signaling pathways. According to the p-value, the top five items from the three categories were selected to plot a histogram (Figure 7B). KEGG enrichment analysis (Figure 7A) indicated that prognostic genes were significantly enriched with arginine and proline metabolism and protein digestion and absorption, etc.

Gene set enrichment analysis was applied to determine their related signaling pathways (Figures 7C–G). COL10A1 was significantly enriched with olfactory transduction and nitrogen metabolism pathways, etc. CTHRC1 was significantly enriched in olfactory transduction and metabolism of xenobiotics by cytochrome p450 pathways, etc. MFAP2 was significantly enriched with olfactory transduction and fatty acid metabolism pathways, etc. P4HA3 was significantly enriched with ribosome and nitrogen metabolism pathways, etc. FAP was significantly enriched in nitrogen metabolism and metabolism of xenobiotics by cytochrome p450 pathways, etc. The mutation site and structure of the five core genes were shown in Figure 8.

Differential expression of the prognostic genes between normal and STAD-related tissues was verified. Results demonstrated that COL10A1, MFAP2, CTHRC1, P4HA3, and FAP were significantly upregulated in STAD-related tissues compared with normal tissues (Figure 9A). Additionally, immunohistochemistry staining of five core genes in STAD and normal tissues was acquired from the Human Protein Atlas database, which showed that differential expression of protein was consistent with gene expression (Figure 9B). However, the immunohistochemical images of COL10A1 were not found. Then, a nomogram (Figure 10A, including TNM stage, age, gender, histologic grade, and risk score) was created to predict the survival rate of patients with STAD at 1, 3, and 5 years. It was found that high total points predicted low 1-, 3-, and 5-year survival rates; however, a low total points did the opposite. The nomogram calibration plots (Figure 10B) indicate that the nomogram was well-calibrated, with mean predicted probabilities for 1- and 3-year OS close to observed probabilities.