STAD-related gene expression data and full clinical annotation with public access were searched in Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) and The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/) database. Patients for whom survival information was lacking were omitted from further evaluation. Six online eligible STAD cohorts, including five GEO STAD cohorts (GSE15459, GSE62254, GSE84437, GSE38749, and GSE34942) and one TCGA-STAD, were selected in this study for further analyses. We excluded data lacking corresponding survival information and negative survival days. Additionally, information was also collected for a cohort from the Affiliated Hospital of Youjiang Medical University for Nationalities (AHYMUN). The AHYMUN (Guangxi, China) cohort included 104 patients diagnosed with STAD at the Department of Hepatology, Affiliated Hospital of Youjiang Medical University for Nationalities, between June 2009 and August 2018. Clinicopathological data were collected from pathology reports or electronic medical records. Samples of STAD and normal gastric tissues were collected during surgery and then processed and stored at the AHYMUN Tissue Bank before the experiments. The research scheme was approved by the research ethics committee of Youjiang Medical University. The raw CEL files were downloaded and a robust multiarray averaging method was adopted with the affy and simpleaffy packages to enable background adjustment and quantile normalization. For the TCGA datasets, RNA sequencing data (FPKM value) of gene expression were downloaded from the Genomic Data Commons (GDC; https://portal.gdc.cancer.gov/) and transformed into transcripts per kilobase million (TPM) values. Three datasets were merged and batch effects from non-biological technical biases were removed using the ComBat algorithm.

A total of 73 HP-related genes were obtained from the GSEA-MSigDB website (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp), where the keywords were Helicobacter pylori infection, including HP_HELICOBACTER_PYLORI_INFECTION and KEGG_EPITHELIAL_CELL_SIGNALING_IN_HELICO:ELICOBACTER_PYLORI_INFECTION. After HP infection, these 73 genes changed. Accordingly, we explored the possible potential biomarkers. The full details of these genes are shown in Table S1 . Non-negative matrix factorization (NMF) analysis was employed to identify distinct HPClusters based on HP-related gene enrichment. R package was utilized to perform an NMF analysis of the HP-related genes. To explore the distinct characteristics in HP-related genes in biological processes, gene set variation analysis (GSVA) was conducted with the hallmark gene set (c2.cp.kegg.v7.2) derived from the MSigDB database. The clusterProfiler package was utilized to perform GO enrichment analysis for the HP-related genes to explore the enrichment features of these genes in different biological processes, cellular components, and molecular functions.

The relationships between HPClusters and clinicopathological characteristics, recurrence-free survival, and prognosis were investigated to examine the clinical value of the two molecular subtypes. Characteristics such as age, sex, race, residual tumor, TNM stage, tumor stage, tumor grade, primary therapy outcome, and future status, were included in the analyses. The overall survival (OS) of the two HPClusters was assessed using Kaplan-Meier curves generated by the “survival” and “survminer” R packages. The clusterProfiler package was used to perform Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment, and Disease Ontology (DO) analyses.

The empirical Bayesian approach of “limma” R package with a fold-change of 1.5 and an adjusted value of p < 0.05 was employed to identify the DEGs between the different HPClusters. A total of 706 DEGs were achieved ( Table S2 ).

A novel risk signature score was constructed to quantify the genetic characteristics of an individual tumor. To build the signature, the DEGs were identified from the HPClusters for LASSO analysis and then multivariate Cox regression analysis. Those DEGs determined as having a significant association with the STAD were selected. Second, patients with STAD were classified into distinct gene clusters using NMF analysis for further analyses. Finally, samples from the TCGA, GSE38749, and GSE84437 cohorts were randomly grouped into the training (n = 399) and testing (n = 399) sets at a ratio of 1:1 and the training set was used to construct the risk signature. Based on HP-related genes, the LASSO and Cox regression algorithm was used to minimize the risk of overfitting using the “glmnet” R package. The change trajectory of each independent variable was analyzed and a 10-fold cross-validation was utilized to construct the signature. Candidate genes were selected using multivariate Cox analysis for risk signature construction in the training set.

The risk signature was defined as follows:

where Expi and Coefi denote the expression and risk coefficient of each gene, respectively. According to the median score, patients in the training set, testing set, and all sets were divided into the low- and high-risk score groups. A Kaplan-Meier survival analysis was then conducted of each group and receiver operating characteristic (ROC) curves were generated.

CIBERSORT was employed to quantify the enrichment of infiltrating immune cells in the TME in the high- and low-risk score groups. The correlation between the immune cells and three selected genes was explored. Moreover, the association between the two groups and microsatellite instability (MSI) was also investigated. An Estimated STromal and Immune cells in MAlignant Tumour tissues using Expression data (ESTIMATE) score was also calculated to examine the purity of the tumor.

Subsequently, the tumor mutation burden (TMB) score was calculated for each STAD patient to investigate the somatic mutations between the high- and low-risk score groups. The differences in expression level of immune checkpoints between the high- and low-risk score groups were visualized using boxplots to predict the response to immunotherapy.

A predictive nomogram was developed using the “rms” package based on the independent prognosis analysis outcomes. In the nomogram scoring system, each variable was matched with a score, with the total score obtained by adding the scores across all variables of each sample. Next, ROC curves for sensitivity and specificity of distinct clinical characteristics and risk scores of 1-, 3-, and 5-year survival events were generated to identify the nomogram with the highest sensitivity.

To test the stability of the signature, we substituted it into the GSE15459, GSE62254, GSE84437, GSE38749, GSE34942, and TCGA cohorts. The sample sizes of the GSE38749 and GSE34942 cohorts were too small (<30) so that the deviation would be too large. Therefore, the GSE84437, GSE38749, and GSE34942 cohorts were mixed for analysis. Furthermore, we searched PubMed for four signatures for STAD within a year (Shao et al., 2021; Wang and Liu, 2021; Zhang et al., 2021; Xiao et al., 2022) and compared them with our signature. Finally, immunohistochemical experiments were performed on the three genes in the obtained signature, and the survivals of the above three genes were verified via the AHYMUN cohort.

All statistical analyses were performed using R version 4.1.0. Statistics with p values < 0.05 were considered statistically significant.

Among the 73 HP-related genes included in this study, a summary analysis of somatic mutations in the TCGA-STAD cohort demonstrated a considerably high mutation incidence of HP-related genes under HP infection conditions ( Figure 2A ). In the TCGA cohort, 151 of 433 samples (34.87%) had noticeable HP-related gene alterations. Among the mutated genes, PTPR71 had the highest mutation frequency, followed by TJP1, PLCG1, EGFR, ATP6V1B1, and PLCG2, while no alteration was observed in CARMIL2, IKBKG, CXCL1, or others.

Next, the copy number variations (CNV) in these HP-related genes were investigated, which indicated a prevalent alteration ( Figure 2B ). Among them, CNV increases were the most prominent in genes including GIT1, ATP6V1C1, CCL5, PAK1, and ATP6V0A1, while CNV decreases were the most significant in ATP6V0E2, CASP3, ATP6V0A4, PTPRZ1, and CDC42. The distributions of the CNV alterations in the 73 HP-related genes on human chromosomes are demonstrated in Figure 2C . To explore the interaction between the above genes and their roles in STAD development and progression, an HP-related gene interaction network was generated ( Figure 2D ). Most HP-related genes (n = 44) functioned as a favorable factor toward STAD, while ATP6V0D2, ADAM10, FLT1, and another 29 genes were considered the risk factors.

In the combination of the GSE38749, GSE84437, and TCGA-STAD cohorts, an NMF analysis was conducted to classify the STAD samples into different molecular subtypes based on the expression of 73 HP-related genes. Two clusters were eventually identified, including 456 samples in the first subtype and 342 in the second ( Figures 3A, B ; Table S2 ). Unsupervised clustering of 73 HP-related genes in the combined cohort also revealed two distinct HP-related modification patterns with significant molecular and clinical characteristics ( Figure 3C ). We termed these two clusters HPClusters C1 and C2, respectively. A prognostic analysis of the two HPClusters indicated a significant survival advantage of the HPClusters C2 ( Figure 3D ). Markedly distinct KEGG pathway enrichment characteristics were also noticed between the two clusters in the GSVA enrichment analysis ( Figure 3E ). In addition, the two HPClusters also showed prominently distinct ICICs ( Figure 3F ). A prevalent amplification in immune cells, such as activated B cells, activated CD8 T cells, macrophages, and natural killer cells, was exhibited by the C1 cluster suggestive of inflammation-promoting characteristics, while activated CD4 T cells, neutrophils, and type 17 T helper cells were enriched in cluster C2. Interestingly, major histocompatibility complex class 1 (MHC-1) expression was significantly higher in the C2 cluster than in the C1 cluster, which is a possible explanation for the higher survival rate in the former than in the latter.

To further understand the characteristics and bio-behaviors of each HPCluster, the clusterProfiler package was employed to perform a GO enrichment analysis of the 73 HP-related genes ( Figure 4A ). These genes showed significant enrichment of biological processes related to cell–substrate adhesion and regulation of the cellular response to growth factor stimuli, which could partly explain the high malignancy and recurrence rates of STAD. Interestingly, these HP-related genes were also significantly enriched in the muscular system processes and muscle contractions. To explore the correlations between the HP-related genes and other diseases and malignancies, a DO analysis was performed, and the result showed that these genes were closely associated with musculoskeletal system cancer, consistent with the above findings ( Figure 4B ). Enrichment of the biological pathways was investigated using the KEGG analysis, demonstrating prominent amplification of pathways related to cell growth, cell proliferation, and angiogenesis ( Figure 4C ). The pathways related to multiple diseases, especially cardiomyopathy, were enriched in the C1 cluster, while those related to wound repair were amplified in the C2 cluster, corresponding with the previous findings.

To further investigate the genetic characteristics and potential biological function of each HPCluster, 706 DEGs were determined between the two HPClusters using the limma R package ( Table S3 ). An NMF analysis was conducted on the 706 DEGs, and two completely distinct gene clusters were identified and termed A and B ( Figures 5A, B ; Table S4 ).

Considerably distinct clinical features are visible in the heatmap in Figure 5C . A significant survival advantage in gene cluster A was demonstrated in the subsequent survival analysis ( Figure 5D ). We compared the expression of HP-related genes in the two gene clusters ( Figure 5E ). Marked differences were found; most focused on higher expressions of favorable factors for STAD in gene cluster A, such as ADAM10, ATP6V0B, and ATP6V1C1, while the expression of risk factors such as ATP6V0E1 and ATP6V1E1 were higher in gene cluster B. This finding corresponds with the former discovery of a better prognosis in gene cluster A.

LASSO and multivariate Cox regression analyses were conducted and three core genes (CRTAC1, BATF2, and CTHRC1) were selected from the 706 DEGs to establish the risk signature ( Figures 6A, B ; Table S5 ). Caret R package was employed to randomly categorize the total STAD patients (n = 798) into two identical training (n = 399) and testing (n = 399) groups, in which the former was used to construct the signature.

According to the results of the multivariate Cox regression analysis, the risk score was constructed as follows:

A significant difference in risk scores was revealed between the distinct gene clusters ( Table S6 ). Patients with a risk score lower than the average were categorized into the low-risk group (n = 177), while those with a risk score higher than the average were grouped into the high-risk group (n = 165). The scoring system was applied to all STAD samples, and patients classified into gene cluster B and HPCluster C1 had markedly higher risk scores than those in gene cluster A and HPCluster C2 ( Figures 6C, D ). An alluvial graph was generated to visualize the distribution of patients in HPClusters C1 and C2, gene clusters A and B, risk score, and future status, from which a higher survival possibility was noticed in the low-risk group ( Figure 6E ). To further investigate the relationships between risk score, molecular and genetic classification, prognosis, and clinical features, a heat map was generated for illustration ( Figure 6F ). This finding verified the above results that gene cluster A and HPCluster C2 had a better prognosis than gene cluster B and HPCluster C1. To better understand the relationship between the risk score and genetic behaviors, the expression of HP-related genes in the low- and high-risk groups were compared and visualized ( Figure 6G ).

To ensure convincing validation of the prognostic capability of the risk signature, a survival analysis was performed in the training and testing groups as well as the originally merged cohort and a similar result of survival advantages in the low-risk group were observed among analyses conducted in all sets ( Figures 7A–C ). In further analyses, favorable HP-related gene expression was higher in the low-risk group, while the opposite result was noted in the high-risk group. It was assumed that a lower risk score was associated with a higher survival possibility and that the survival rate would decrease as the risk score increased; this hypothesis was subsequently verified by the distribution plot of the risk score ( Figures 7D–I ). To assist in the comprehension of the influence of the three hub genes on oncogenesis, how the three core genes were distributed in the low- and high-risk groups within the training group, the testing group and total samples were illustrated via heatmaps ( Figures 7J–L ). To further verify the accuracy and efficiency of the risk signature, ROC curves in the training, testing, and combined groups were generated ( Figures 7M–O ). Analysis of the 1-, 3-, and 5-year prognostic categorization and prediction efficiencies demonstrated considerably higher AUC values of the risk signature, which substantiated the ability of the risk signature to predict the future status of STAD patients.

In addition, a survival analysis was conducted in groups of different clinical features such as age (>60 years, <60 years) and sex with no results challenging the higher survival rate of the low-risk group. The above results suggest that lower risk scores were associated with higher survival possibilities ( Figures 8A–D ). Subsequent univariate and multivariate Cox regression analyses of the merged cohort also confirmed that the risk signature could act as an independent prognostic model for STAD ( Figures 8E, F ). To testify to the ability to predict future status, risk scores of STAD samples from alive versus dead patients were compared and visualized, and the former showed significantly lower risk scores ( Figure 8G ). Analyses of the future statuses of the low- and high-risk score groups revealed a prominently higher survival rate in the low- versus high-risk group (63% vs 46%, respectively) ( Figure 8H ).

CIBERSORT, Tumor Immune Estimation Resource, and other algorithms were employed to evaluate the association between risk score and immune cell infiltration, which was visualized via scatter diagrams and a heat map ( Figure 9A ). The analysis demonstrated that the risk score was positively correlated with macrophage M2, resting mast cells, monocytes, and resting memory CD4 T cells. Other immune cell types, such as M1 macrophages, resting NK cells, activated memory CD4 T cells, CD8 T cells, and follicular T helper cells showed negative correlations. However, no significant difference was seen between the low- and high-risk groups despite prominently lower stromal and ESTIMATE scores in the low-risk group ( Figures 9B, C ). Data of the STAD samples were analyzed further using the Tumor Immune Dysfunction and Exclusion algorithm. STAD patients in the low-risk score group showed higher MSI expression and significantly obvious T cell infiltration, indicative of better outcomes of ICB immunotherapy ( Figure 9D ). In the ESTIMATE analysis, significant differences were detected in stromal and ESTIMATE scores between the low- and high-risk score groups ( Figure 9E ). The high-risk score group exhibited prominently higher existence of tumor cells and higher tumor purity than the low-risk score group. In addition, the risk score and cancer stem cell (CSC) index values were synthesized to evaluate the correlation between the risk signature and CSC in STAD ( Figure 9F ). It was revealed that a negative linear correlation existed between the risk score and CSC index (R = -0.59, p < 2.2^-16), indicating that the lower the risk score, the more distinct the stem cell properties and lower the cell differentiation degree. Subsequently, patients with lower risk scores showed obviously higher expressions of PD-1, PD-L1, and CTLA-4, suggesting a potential beneficial response in the low-risk score group to ICB therapy ( Figures 9G–I ).

Further analysis investigated the relationship between immune infiltration and the three HP-related core genes. As Figure 9E shows, most immune cell types, especially activated memory CD4 T cells, showed a significantly positive correlation with BATF2, while the slight opposite occurred with CRTAC1 and CTHRC1. The distribution differences of somatic mutation between the low- and high-risk score groups were then analyzed using the maftools package ( Figures 10A, B ). To further understand how TMB and risk score can predict the prognosis of STAD patients, more survival analyses were performed on groups with different TMB and risk scores. The low-risk score group with high TMB had the highest survival probability ( Figures 10C, D ). The correlations between TMB and risk score and between the HPClusters and gene clusters were then evaluated ( Figures 10E, F ). This result further verified that a low-risk score was linked to a high TMB. Clusters C2 and B also showed significantly higher TMB than the other two clusters. Consistent with the above findings, a more extensive tumor mutation burden was seen in the low-risk score group than in the high-risk score group, with overall mutation rates of 89.83% versus 86.06%, respectively, which was subsequently confirmed by the TMB quantification analyses ( Figure 10G ). To further assess the ability of the risk signature to predict STAD patients’ response to ICB therapy, we analyzed the MSI in the low- and high-risk score groups ( Figure 10H ). The low-risk score group contained a higher percentage of high MSI than the high-risk score group, suggesting that the former was more sensitive to immunotherapy and had greater clinical advantages. This correlation between the risk signature and MSI was confirmed further ( Figure 10I ).

A nomogram containing the risk score and clinicopathological parameters was established to predict the 1-, 3-, and 5-year survival rates of STAD patients ( Figure 11A ). A subsequent calibration plot also suggested that the proposed nomogram had a similar performance compared to an ideal model ( Figure 11B ). To perform a more accurate and convincing validation, clinical ROC curves for 1-, 3-, and 5-year survival of STAD patients were generated and the nomogram showed the highest AUC values, suggesting its superior survival predictive ability compared to other parameters ( Figures 11C–E ).

To verify signature reliability, we applied the model formula to cohorts GSE15459, GSE62254, GSE84437+GSE38749+GSE34942, and TCGA-STAD ( Figures 12A–E ). We compared the results with the signatures reported by Zhang et al., Shao et al., Wang et al., and Xiao et al. ( Figures 12F, G ). Our signature had significant effects on the different online cohorts. Moreover, the C-index value of our signature was the highest among the five signatures, and its RMS value was <2, representing low dispersion and high reliability.

The three genes in the signature, CRTARC1, BATF2, and CTHRC1, were analyzed by immunohistochemistry in the STAD samples ( Figure 13A ) and combined with the sample information of the 104 AHYMUN cohorts ( Figure 13B ). As previously analyzed, the high expression of CRTARC1 and CTHRC1 in the STAD samples represented lower OS, while higher BATF2 expression predicted a better prognosis.