The transcriptome data and clinical materials of GC (n=378) were downloaded from TCGA (https://portal.gdc.cancer.gov/). The clinical features are detailed in Supplementary Table S1 . Moreover, we selected three independent validation datasets (GSE15459, GSE62254, GSE84437) from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) and obtained their normalized microarray gene expression data and clinical data. We obtained a list of 300 PRGs from previous literature ( Supplementary Data Sheet 1 ) (20).

After collection and preprocessing the data of GC, the Univariate Cox regression analysis was performed on the PRGs collected to identify PRGs with prognostic value (P<0.05). Genetic mutations of the prognostic PRGs were analyzed on the cBioportal online tool (https://www.cbioportal.org/) using the Stomach Adenocarcinoma (TCGA, Firehouse Legacy) dataset. Moreover, a PPI network diagram of the prognostic PRGs was constructed with the STRING database (http://string-db.org/) and graphed with the Cytoscape software (21) (version 3.7.2). Differential expression analysis for the prognostic PRGs between tumor and normal tissues was performed using the “limma” package in R (version 4.2.3) (22).

The least absolute shrinkage and selection operator (LASSO) regression analysis of the candidate prognostic PRGs was performed to construct a prognostic gene signature. Then, we calculated each patient’s risk score. The calculation formula is as follows:

Based on the median value of the risk score, patients in the TCGA training group were divided into high- and low-risk groups. To further verify the predictive ability of the model, three independent validation datasets (GSE15459, GSE62254, and GSE84437) were included in our study. Risk scores were calculated separately for each sample in the training cohort and the GEO validation cohorts based on the same risk formula. Based on the median risk score, we could divide the patients into two subgroups of high risk and low risk to explore the prognostic differences between the two groups. The Kaplan-Meier curves and receiver operating characteristic curves (ROCs) were constructed for the training cohort and validation cohorts.

To determine if the PLT signature may serve as a standalone predictive factor in patients with GC, we preformed multivariate Cox regression analysis. A nomogram for predicting overall survival (OS) at 1, 2, and 3 years in clinical patients was constructed using the “rms” R package based on the patient’s age, histologic grade, gender, stage and risk scores.

We utilized the “limma” R package (22) to identify differentially expressed genes (DEGs) between the high-risk group and the low-risk group with the criteria of fold change (FC) > 2 and false discovery rate (FDR) < 0.05. To further investigate the function of the DEGs, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were analyzed using hypergeometric distribution testing by the “ClusterProfiler” R package (23). “circlize” R package (24) visualizes the GO and KEGG results. Finally, Gene Set Enrichment Analysis (GSEA) with the Kolmogorov–Smirnov (KS) test was performed to find enriched KEGG pathways, the ridge plot was used to present the details of GSEA via the “ggstatsplot” R package.

The Immuno-Oncology Biological Research (IOBR) R package (25) (version 0.99.9) was used to analyze the immune features and immune cell infiltration in high- and low-risk groups. Based on the 186 TME-associated signatures in the R packet IOBR, we calculated the sample enrichment score. We assessed the expression of common immune checkpoint genes between high-risk and low-risk groups. The relationship between risk score and immune checkpoint genes was analyzed with Pearson correlation test. The CIBERSORT algorithm in the IOBR package was used for calculating the relative abundance of 22 kinds of immune cells in TCGA-GC cohort, and the ESTIMATE algorithm for calculating each sample’s matrix and immune scores.

The Cancer Therapeutics Response Portal (CTRP) database contains data on the sensitivity of different tumor cells to different chemotherapeutic drugs. We employed the R package “oncoPredict” (26) to calculate the sensitivity of individual GC patient to different chemotherapeutic drugs based on the gene expression data (log2(TPM + 1)). Then, the difference in the area under the dose–response curve (AUC) values between high-risk and low-risk groups was evaluated.

The protein expression of ten signature genes in GC and normal tissues was determined using immunohistochemistry (IHC) from the Human Protein Atlas (HPA) (https://www.proteinatlas.org/), which is a valuable database providing extensive transcriptome and proteomic data for specific human tissues and cells.

The normal human gastric epithelial cell line GES-1, and four human gastric cancer cell lines, MKN45, N87, HGC27, and KATO-3, were authenticated by STR profiling. All cell lines were cultured in RPMI-1640 containing 1% penicillin/streptomycin and 10% fetal bovine serum. Cells were grown in 5% CO2 at 37°C. Total RNA was extracted using TRIzol (TransGen Biotech, China). Complementary DNA (cDNA) was synthesized using the GoScript™ Reverse Transcription Mix and Oligo(dT) kit (Promega, United States). Real-time PCR was performed using SYBR Green PCR Master Mix (FastStart Universal SYBR Green Master, Roche). Relative gene expression levels were normalized to the levels of GAPDH using the ΔCt method. Each experiment was operated in technical triplicate. The amplification primer sequences of each gene are detailed in Supplementary Data Sheet 6 .

R software version 4.2.3 was used to conduct the statistical analysis, and p-values and FDR q-values below 0.05 were regarded as statistically significant.

The primary design of this study was depicted in the flow chart ( Figure 1 ). A total of 30 PRGs were significantly associated with prognosis of GC patients based on the Univariate Cox regression analysis ( Supplementary Data Sheet 2 ). Genetic mutations of the 30 prognostic PRGs were analyzed through cBioPortal online tool for GC patients. Genes with a mutation rate no less than 5% are shown in Figure 2A . COL1A2 had the highest mutation rate (13%) among 478 patients, followed by F5 (7%), GNG11 (7%), FN1 (7%), and DGKG (7%). We constructed a protein interaction network for the 30 genes based on the STRING database. As shown in Figure 2B , the FN1, ALB, SERPINE1, COL1A1, COL1A2, and ITGB1 genes were at the core of the protein network interaction. The differential expression analysis of the candidate prognostic PRGs revealed 23 differentially expressed genes, including 18 upregulated and 5 down-regulated genes ( Figure 2C ).

The forest plot of the 30 prognostic PRGs obtained by univariate Cox regression analysis was shown in Figure 3A . Then we constructed a predictive prognostic model consisting of 10 PRGs by LASSO regression analysis ( Figures 3B–D ). They were SERPINE1, ANXA5, DGKQ, PTPN6, F5, DGKB, PCDH7, GNG11, APOA1, and TF. The coefficient and HR value of multivariate Cox regression analysis is shown in the form of a forest map ( Figure 3B ). A linear prediction model was developed based on the weighted regression coefficients of 10 prognostic PRGs, calculated as r i s k   s c o r e = ( − 01357 × S E R P I N E 1   exp ) + ( 0.12273 × A N X A 5   exp ) + ( − 0.0927 ) × D G K Q   exp + ( − 0.0722 × P T P N 6   exp ) + ( 0.05354 × F 5   exp ) + ( 0.04579 × D G K B   exp ) + ( 0.04079 × P C D H 7   exp ) + ( 0.03418 × G N G 11   exp ) + ( 0.03202 × A P O A 1   exp ) + ( 0.02374 × T F   exp ) . Of these, SERPINE1, ANXA5, F5, DGKB, PCDH7, GNG11, APOA1, and TF showed significant positive correlations with risk scores, while DGKQ and PTPN6 showed significant negative correlation with risk scores.

After establishing the predictive prognostic model based on 10 prognostic PRGs for GC, we computed the risk score for each GC patient based on the LASSO coefficients and expression value for each PRG ( Supplementary Data Sheet 3 ). We contrasted the distribution of risk score, the survival status and the heatmap of GC patients in the TCGA cohort ( Figure 4A ). The risk curves and scatter plots implied that mortality was positively related to the risk score in the TCGA cohort. The heatmap unveiled that higher DGKQ and PTPN6 expression were detected in the low-risk group, while the other eight genes (F5, APOA1, TF, ANXA5, SERPINE1, PCDH7, DGKB, and GNG11) were highly expressed in the high-risk group. Kaplan-Meier analysis was used to analyze the survival and prognosis of GC patients in TCGA. As shown in the Figure 4B , patients in the low-risk group had a better prognosis, while patients in the high-risk group had a worse prognosis (P<0.001). The AUCs of 1-year, 2-year, and 3-year survival ROC curves predicted by the PLT signature were 0.670, 0.695, and 0.707, respectively, suggesting the efficiency of PLT signature in predicting prognosis for GC to a certain extent ( Figure 4C ).

To further demonstrate the stability and reliable generalization of our model, the GSE15459, GSE62254, and GSE84437 cohorts were used as the external validation cohorts. The Kaplan-Meier curves showed a significant difference in prognosis between the high-risk and low-risk patients in these three cohorts, respectively, with a more significant survival advantage for patients in the low-risk group (P = 0.001, P = 0.003, P < 0.001, respectively) ( Figures 4D–F ). The ROC curve was used as a tool to predict the survival time of patients at 1-, 2-, and 3- years. The AUCs at 1-, 2-, and 3- years for the GSE15459 cohort were 0.670, 0.633, and 0.662, respectively ( Figure 4G ). The AUCs for the GSE62254 cohort were 0.667, 0.606, and 0.608, respectively ( Figure 4H ). The AUCs for the GSE84437 cohort were 0.599, 0.608, and 0.611, respectively ( Figure 4I ). This indicates that the model has an excellent predictive effect.

To validate the reliability and clinical value of the biological signature constructed based on PRGs as a predictor of prognosis, we conducted multivariate Cox regression analysis including common clinical characteristics ( Supplementary Data Sheet 4 ). It is shown that in the multifactorial cox analysis, tumor stage (P<0.001) and risk score (P<0.001) were all independent prognostic factors significantly associated with patient prognosis ( Figure 5A ). Based on the above analysis, in order to be able to predict patients’ prognosis quantitatively and to inform clinical decision-making, we integrated the risk score and clinical indicators to construct Nomogram plots as a means of predicting the probability of prognostic survival at 1, 2, and 3 years ( Figure 5B ). We then used time-dependent ROC curve analysis to compare the predictive accuracy between the nomogram, risk score, and common clinicopathological features ( Figure 5C ). The results showed that risk score had a much greater AUC value than the rest of the individual clinicopathological features, and the nomogram model suggested higher prognostic accuracy at 1-, 2-, and 3-year OS with a larger AUC than risk score. The time-dependent AUCs of the nomogram for predicting 1-, 2-, and 3-year OS were 0.708, 0.763, and 0.742, respectively. Combined with these results, this suggests that our PLT signature is more practical and influential for clinical decision making and is more suitable as a clinical decision tool for predicting the prognosis of patients with GC in the clinical setting.

We performed DEGs analysis between high-risk and low-risk groups on the TCGA cohort, and the results showed that 2,442 DEGs were differentially expressed between the high-risk group and the low-risk group based on the criteria of P < 0.05. Among that, 2,249 genes were up-regulated, and 193 genes were down-regulated. The volcano plot of DEGs were displayed in Figure 6A . All of the upregulated and downregulated genes were demonstrated in Supplementary Data Sheet 5 . The results of GO analysis can be divided into three categories: biological process, cellular component, and molecular function. Where in biological processes, such as axonogenesis, extracellular matrix organization, extracellular structure organization; Cellular components, such as collagen−containing extracellular matrix and synaptic membrane; And molecular functions, such as extracellular matrix structural constituent, G protein−coupled peptide receptor activity, peptide receptor activity, and glycosaminoglycan binding were significantly enriched ( Figure 6B ). KEGG pathways were enriched in Neuroactive ligand−receptor interaction, cyclic adenosine monophosphate (cAMP) signaling pathway, Calcium signaling pathway, Cell adhesion molecules, and extracellular matrix (ECM)−receptor interaction ( Figure 6C ). Then, the GSEA method was applied to identify the significantly enriched KEGG pathways in the high-risk samples. The ridgeplot showed that several pathways, such as calcium signaling pathway, cAMP signaling pathway, ECM-receptor interaction, focal adhesion, and neuroactive ligand-receptor interaction, were significantly enriched in the high-risk group ( Figure 6D ). DNA replication, base excision repair, homologous recombination, nucleotide excision repair were the pathways that were substantially enriched in the low-risk group. GSEA plot of important pathways enriched in the high-risk group was shown in Figure 6E .

To further elucidate differences in the immune microenvironment of patients between high-risk and low-risk groups, we compared the enrichment scores of TME cells-related signatures between two groups. The results showed that T cell-related signatures [T cell accumulation, T cell exhaustion, T cell regulatory (27)] and tumor-associated macrophage-related signatures [Macrophages Bindea et al. (28), TAM_Peng_et_al (27)] had significantly higher enrichment scores in the high-risk group compared to the low-risk group ( Supplementary Figure S1 ). Additional examination of TME signatures employing the IOBR package unveiled an immunosuppressive, exclusive, and exhausted TME in the high-risk group ( Supplementary Figures S2A–C ). Furthermore, patients in the low-risk group exhibited higher scores in DNA damage response (DDR), mismatch repair, and homologous recombination ( Supplementary Figure S2D ), suggesting that they may be more sensitive to immunotherapy. High-risk patients demonstrated more pronounced epithelial-mesenchymal transition (EMT) signatures ( Supplementary Figure S2E ). Taken together, these findings suggest an immunosuppressive TME in the high-risk group. The extent of immune cell infiltration in patients in the TCGA cohort was then assessed. The results of ESTIMATE suggested that stromal score, and ESTIMATE score were higher in the high-risk group ( Figure 7A ). We then estimated the proportion of 22 types of immune cells in each sample by CIBERSORT algorithm. The difference in the proportion of each type of immune cell between two risk groups was shown Figure 7B . The results revealed that compared with the low-risk group, memory B cells (P < 0.001), follicular helper T cells (P < 0.001) exhibited lower infiltrating levels in the high-risk group. However, samples in the high-risk group had a significant increase in the fraction of naïve B cells (P < 0.01), monocytes (P < 0.001) and macrophages M2 (P < 0.01). We also explored the relationship between risk score and common immune checkpoint genes, including programmed cell death 1 (PDCD1), PDCD1 ligand 1 (PDCD1L1/CD274), cytotoxic T-lymphocyte-associated antigen 4 (CTLA4), PDCD1 ligand 2 (PDCD1LG2), hepatitis A virus cellular receptor 2/T-cell immunoglobulin mucin receptor 3 (HAVCR2/TIM3), lymphocyte activating 3 (LAG3), and T cell immunoreceptor with immunoglobulin and ITIM domain (TIGIT). As displayed in Supplementary Figure S3A , the levels of HAVCR2 and PDCD1LG2 were increased in the high-risk group in GC. However, no significant correlation was observed between the immune checkpoint molecules and the risk score ( Supplementary Figure S3B ).

In order to find more effective drugs for patients in the high-risk group, we further studied the sensitivity of tumor cells to chemotherapeutic drugs between different risk groups based on CTRP database. The AUC value represents the degree of drug sensitivity. An increasing AUC value represents a lower drug sensitivity. We found that patients in the high-risk group tended to be less sensitive to oxaliplatin, doxorubicin, and mitomycin, but more sensitive to thalidomide, MK-0752, and BRD-K17060750 ( Figure 8 ).

IHC results of the protein expression of the signature genes from HPA database were displayed in Figure 9 . SERPINE1 was expressed at a low level in stomach normal tissues and was not detected or expressed at a low level in tumor tissues. ANXA5 was expressed at a low to medium level in both tumor and stomach normal tissues. DGKQ was not detected or expressed at a low level in GC tumor tissues but had medium to high expression levels in stomach normal tissues. PCDH7 was expressed at a high level in stomach normal tissues and had diverse expression levels in GC tumor tissues, ranging from low expression and medium to high expression. GNG11 was not detected in stomach normal tissues but had various expression levels in GC tumor tissues, which was from not detected and low to medium expression. APOA1 was not detected in both tumor and stomach normal tissues. TF was not detected or expressed at a low level in normal tissues and had diverse expression levels in GC tumor tissues, which was from not detected and low expression to medium expression. Data for F5 and DGKB were lacking and therefore not presented.

The qRT-PCR analysis was performed to further verify ten signature genes in normal and tumor cells ( Figure 10 ). The results showed that the expression of PTPN6, F5, DGKB, PCDH7, and TF was elevated in GC cell lines, whereas the expression of SERPINE, ANXA5, DGKQ, and GNG11 was decreased. APOA1 expression was not detected in the GES-1 and four human gastric cancer cell lines.