The clinicopathological information and corresponding gene expression data of GC patients were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). A total of 684 cases (GSE13861, GSE29272, GSE62254, and GSE26942) were examined in the entire cohort. 330 TCGA-STAD samples and their corresponding clinicopathological data were extracted from the TCGA (https://portal.gdc.cancer.gov/) to be used for the external validation cohort. A list of focal adhesion-related genes was retrieved from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb) to aid this analysis.

Gene symbols for each gene matrix file were extracted based on the corresponding platform file used by the Perl software. The batch effect was adjusted using the Empirical Bayes method (“sva” package) among the series. Finally, the entire cohort was randomly divided into a training and a test cohort in a 7:3 ratio.

A total of 199 focal adhesion-related genes were selected from the MSigDB database. The relationship between these genes and gastric cancer prognosis was investigated through univariate Cox, LASSO regression, and multivariate Cox analyses using R packages “glmnet” and “survival.”

We calculated the risk score of every patient using the regression coefficient value and expression of each gene, as follows: = Σβi × Expi, where Expi refers to the gene expression level of the focal adhesion-related genes, and β is the LASSO Cox regression coefficient of the corresponding gene. Based on their risk scores, samples were categorized as high- or low-risk. The prediction accuracy of the signature was evaluated using time-dependent ROC curves and Kaplan-Meier survival analysis.

Based on the focal adhesion risk score, the patient’s gender, age, and American Joint Committee on Cancer AJCC stage, a predictive nomogram was developed to accurately predict the OS of GC patients. The precision of the nomogram was evaluated using calibration and receiver operating characteristic (ROC) curves.

The GSEA program was obtained from the GSEA website (http://www.gseamsigdb.org/gsea/index.jsp) to identify the functional enrichment pathways regulated by the focal adhesion -related genes signature. Additionally, the “c2.cp.kegg.v7.4.symbols.gmt” and “c5.go.v7.4.symbols.gmt” gene sets were extracted from the molecular signatures database and used as the target enrichment sets for GSEA analysis.

The Estimation R package was used to estimate the immune and stromal component scores in the tumor microenvironment (TME) of each GC sample. The immune and stromal scores were then calculated and displayed, along with an estimated score. A positive correlation between the risk score and the immune, stromal, and overall scores suggests a higher risk score and a greater proportion of corresponding TME components.

The R packages “survival” and “survminer” were used to analyze the TME scores. A total of 684 samples with available survival data were divided into high- and low-scoring groups based on the median values of the immune, stromal, and estimated scores for subsequent survival analysis.

The single sample gene set enrichment analysis (ssGSEA) in the “gsva” R package quantifies the infiltration statuses of 16 immune cells and the activities of 13 immune-related pathways in the high- and low-risk

The drug sensitivity analysis was conducted using cellminer (https://discover.nci.nih.gov/SclcCellMinerCDB/) Database data, screening FDA-approved and clinical trial data, and analyzing the relationship between focal adhesion-related genes expression level and drug sensitivity. Spearman’s correlation analysis was conducted to determine the correlation using R software, and the top 16 drugs were selected.

The accuracy of the results was further validated using the STAD data from the KM-plotter database (https://kmplot.com/).

Download the GSE112301 dataset from the GEO website and use the CreateSeuratObject function from the SeuratR package to create a Seurat object that contains the basic information of the single-cell dataset. Next, quality control is performed on the data, including filtering out low-quality cells, followed by reducing the dimensionality of the data using principal component analysis (PCA). Visualize the PCA results to better distinguish the differences between cells in different tissue. Finally, use the FindConservedMarkers function of Seurat to analyze the gene expression of GC samples and normal samples.

The statistical analyses were conducted using R software version 4.0.0, and KM survival analysis assessed the differential OS durations between the high- and low-risk groups. p < 0.05 was set as the significance threshold.

As illustrated in Figure 1, upon exclusion of cases with a survival time of fewer than 30 days and normal cases, 684 samples were collated in the four GEO datasets (GSE13861, GSE26942, GSE29272, and GSE62254). These cases were randomly divided into a training (478) or test cohort (206) in a 7:3 ratio. Table 1 summarizes the patient clinical characteristics. The “limma” package extracted genes associated with focal adhesion in the GEO database. In the training cohort, the univariate Cox and LASSO regression analyses were conducted to screen eight genes associated with GC patient OS, as depicted in Figure 2A. The risk score was calculated according to a linear combination of the expression levels of the eight focal adhesion-related genes and corresponding regression coefficients (Table 2).

Risk score = CMOP × (0.3294) + FLNC × (0.2545) + ITGB5 × (0.7737) + LAMC1 × (0.8969) + TBHS1 × (0.5076) + THBS4 × (0.0156) + VEGFB × (0.8368) + VWF × (0.2495).

Based on the optimal cut-off value of the risk score determined using the Survminer R package, the patients were categorized into either a high- or low-risk group. Further independent prognostic analysis of all the key genes revealed a significant association between high expression levels of these genes and poor prognosis in the training cohort (Supplementary Figure S1).

The distribution of the risk scores and the survival statuses of patients in the training cohort are displayed in Figures 2B–D. KM analysis determines the differences between the two groups in the training cohort (Figure 3A). Lastly, time-dependent ROC exhibits the prognostic values of our signature (Figure 3D).

Test and external validations were conducted. Consistent with the training cohort results, OS was lower in high-risk patients than in low-risk patients (Figure 3B, p < 0.001 in test cohort; Figure 3C, p = 0.007 in TCGA cohort). The area under the ROC curve demonstrated that the signature could precisely predict GC prognosis (Figures 3E, F). The AUCs for the test and external validation cohorts were 0.667 and 0.627 at 3 years, respectively, and 0.652 and 0.701 at 5 years, respectively. The accuracy of the signature was evaluated using PCA and t-SNE analyses. Both PCA and t-SNE plots revealed that the high- and low-risk groups had different directions in the training (Figures 3G, J), test (Figures 3H, K), and external cohorts (Figures 3I, L).

The Cox regression analysis demonstrates the relationship between the risk scores acquired from the prognostic model with other clinical parameters. Based on the uni- and multivariate regression analyses, stage (p < 0.001, HR = 3.369; p < 0.001, HR = 3.050, respectively) and risk score (p < 0.001, HR = 3.314; p < 0.001, HR = 2.980, respectively) were independent OS prognostic factors in the training cohort (Figures 4A, B). The test cohort was validated, whereby both stage (p < 0.001, HR = 3.248; p < 0.001, HR = 2.903, respectively) and risk score (p < 0.001, HR = 3.025; p = 0.002, HR = 2.676, respectively; Figures 4C, D) were demonstrated to be independent risk factors for OS in GC patients.

To investigate the prognostic value of the developed model in different patient populations based on their clinical characteristics, the training cohort was further divided into subgroups, and OS was estimated between high- and low-risk groups in each subgroup. The KM analysis showed that the risk score could distinguish differences between various subgroups, such as age, gender, and stage (Figures 5A–F). Similarly, in the test cohort, the high- and low-risk groups demonstrated differences in age and gender, but no significant differences were observed in the stage I-II subgroup, which may be attributed to the relatively small sample size in our study (Supplementary Figures S2A–F).

The relationship between risk score and corresponding clinical characteristics of GC patients (n = 684) in the GEO dataset was analyzed. The risk scores of stage I + II patients were significantly lower than those of stage III + IV patients (p = 1.7e-09, Figure 6A). Similar results were shown in the age subgroup (p = 0.0095, Figure 6B). No significant relationship between the risk score and gender (p = 0.72, Figure 6C).

To accurately predict the prognosis of GC patients, a nomogram was developed based on uni- and multivariate regression analyses to predict 1-, 3-, and 5-year OS rates (Figure 7). Additionally, 3 years time-dependent ROC analysis revealed that the sensitivity of the nomogram was higher than other clinicopathological features in training cohort (Figure 8A), test cohort (Figure 8B), entire cohort (Figure 8C) and external cohort (Figure 8D). 5 years time-dependent ROC analysis for nomogram in training cohort (Figure 8E), test cohort (Figure 8F), entire cohort (Figure 8G) and external cohort (Figure 8H). The calibration plots for the training, test, entire cohort and external cohort were in agreement between the actual OS and the predicted from the nomogram (Figures 8I–8L).

The nomogram generated to predict the 1-, 3-, and 5-year OS rates of GC patients was found to be accurate, as evidenced by the 3-year and 5-year time-dependent ROC analysis and the calibration plots. These results suggest that the nomogram is a reliable tool for predicting the prognosis of GC patients.

To explore the functional and signaling pathway differences between the high- and low-risk score groups, GSEA was performed on the gene sets “c5.go.v7.4.symbols.gmt” and “c2.cp.kegg.v7.4.symbols.gmt.” The top five pathways and gene functions in the high- and low-risk groups are displayed in Figures 9A–D.

The tumor microenvironment (TME) plays a critical role in regulating tumor treatment resistance and is associated with tumor occurrence, development, and metastasis. It includes various components such as tumor cells, immune cells, stromal cells, and a variety of cytokines Changes in the TME, including alterations in the immune cell components, can promote tumor progression. To analyze the distribution of immune cells in the GC TME and investigate the interaction between GC tumors and immune cells, we utilized the ssGSEA tool to predict 16 common immune cells and 13 immune-related functional components based on GC gene expression profile data.

The low-risk group of patients showed higher ratios of B cells, T regulators, and follicular helper T cells than the high-risk group (Figure 10A). In addition, the low-risk group exhibited higher levels of antigen-presenting cell (APC) co-inhibition, inflammation promotion, and T cell inhibition compared to those with a high-risk score (Figure 10B). Additionally, we demonstrated that the expression levels of CD200, CD28, CD40, CD44, CD86, LAIR1, NRP1, TNFRSF4, TNFRSF8, TNFSF18, TNFSF4, and VTCN1 in the high-risk group were higher than those in the low-risk group (Figure 10C). The findings suggest that the immune microenvironment may be partly associated with the OS prognosis of GC patients with high expression of focal adhesion-related genes.

N6-methyladenosine (m6A) is the most abundant RNA modification in eukaryotic cells (Yue et al., 2015). Extensive RNA processing and metabolism research revealed that m6A is a key contributor to cancer development. m6A is a potential prognostic marker involved in multiple aspects of cancer treatment (Ma et al., 2019). To assess the relationship between m6A expression and our GC prognostic signature, the levels of 13 m6A genes in different GC samples were estimated. It is found that an elevated expression of FTO, METTL3, YTHDC1, and YTHDF1 genes in the high-versus low-risk group (Figure 10D).

TME consists of diverse immune and stromal cells linked to disease development, prognosis, and treatment outcome. Based on our ESTIMATE algorithm, TME was separated and scored into stromal, immune, and estimate subcomponents to investigate potential relationships between this study’s risk scores and TME. A high immune or matrix score indicates a high proportion of the immune or matrix components in the TME. The ESTIMATE score is the sum of the immune and stromal scores, indicating the combined proportion of these two components in TME. In our study, patients in the high-risk group in entire cohorts had higher stromal, immune, or ESTIMATE scores (Figure 11A) than those in the low-risk group.

To further investigate the impact of different components of the TME on GC patient survival, the entire cohort was divided into subgroups based on the median immune, stromal, and ESTIMATE scores as cutoff points. As shown in Figures 11B, D, patients with high stromal and ESTIMATE scores had worse overall survival than those with low stromal and ESTIMATE scores (p < 0.001, p = 0.01, respectively). However, the survival rate was similar between patients with high and low immune scores (p = 0.933) (Figure 11C). To better understand the association between the immune microenvironment and GC prognosis, a heatmap was generated to display the distribution of immune cell scores in the high- and low-risk groups, as shown in Figure 12.

The correlation between drug Z score and genes was analyzed, and the top 16 significant drug-gene pairs are displayed in Figure 13. A total of 246 drugs showed statistical differences, as shown in Supplementary Table S2. Among them, Dasatinib, XAV-939, and Staurosporine exhibited the most positive correlation with hub gene expression. In contrast, Palbocic, Oxaliplatin, and Ribavirin were negatively correlated with the expression of hub genes (Figure 13).

To investigate the prognosis of patients with different hub genes expressions, the clinical data of STAD in the KM-plotter database were analyzed. Patients with high hub genes expression had better overall survival (OS) than those with low expression, except for THBS1 (Figure 14).

The Single-cell RNA sequencing dataset GSE112302 was used for further analysis in high resolution. A total of three GC and three normal samples were included in our study.

After quality control of the data, standardization and normalization were performed, followed by PCA and UMAP (Supplementary Figures S3A, B). A total of 305 normal cells and 401 tumor cells were included in the analysis. The UMAP plots of each gene in different tissues are displayed in Supplementary Figure S3C. Most genes exhibited high expression in tumor cells, except for ITGB5, which may be due to the limited sample size (Figure 15).