Helicobacter pylori virulence gene expression defines a stage-associated transcriptional signature in gastric cancer

Abstract

BACKGROUND

Helicobacter pylori (H. pylori) is an important risk factor for gastric adenocarcinoma; however, clinical outcomes differ substantially according to the virulence profile of the infecting strain. Diagnostic approaches currently used generally determine pathogen status in a binary manner and therefore do not capture complex molecular variation among strains or variation associated with different disease stages.

AIM

To determine whether H. pylori virulence gene expression varies across clinical stages of gastric cancer (GC) and whether such variation can be used to establish a transcriptional signature associated with disease stage.

METHODS

A two-tier molecular framework was applied to 400 participants recruited in Pakistan, comprising 200 GC cases and 200 non-malignant controls. Tier one examined the prevalence and risk associations of four virulence genes: Cytotoxin-associated gene A (cagA), vacuolating cytotoxin A (vacA), blood group antigen-binding adhesin A (babA), and induced by contact with epithelium gene A (iceA). Tier two quantified transcript expression in 159 H. pylori-positive GC cases spanning tumor-node-metastasis stages I-IV. Relationships between expression and clinical stage were analyzed by ordinal logistic regression, principal component analysis, and a recursive partitioning (decision tree) classifier. Receiver operating characteristic analysis and cross-validation were used to evaluate model performance.

RESULTS

Tier one demonstrated a distinct hierarchy among virulence factors: vacA had the strongest association [odds ratio (OR) = 22.9], with cagA ranking second (OR = 18.7). Tier two identified stage-related differences in virulence gene expression, most notably a decline in babA expression, while cagA and vacA remained comparatively elevated at later disease stages. In multivariate analysis, babA downregulation (OR = 0.536, P < 0.001) and iceA expression (OR = 0.802, P = 0.010) emerged as the most significant independent predictors of clinical stage. In internal validation, the combined four-gene signature yielded an area under the curve (AUC) of 0.769 (95%CI: 0.687-0.851) for overall stage discrimination. This apparent AUC fell to a pooled 5-fold cross-validated AUC of 0.747 (95%CI: 0.664-0.823) on held-out folds, indicating good generalizability. The decision-tree classifier separated early stages (I-II) from advanced stages (III-IV) with an apparent AUC of 0.867 (95%CI: 0.802-0.932; sensitivity 74.5%, 95%CI: 61.6-84.2%; specificity 81.7%, 95%CI: 73.2-88.0%) on the full cohort, which fell to a nested cross-validated AUC of 0.68 (95%CI: 0.59-0.77; sensitivity 54.5%, 95%CI: 41.5-67.0%; specificity 75.0%, 95%CI: 65.9-82.3%) on held-out folds.

CONCLUSION

Expression of H. pylori virulence genes differs across stages of GC. This hypothesis-generating transcriptional signature demonstrated superior discriminatory performance compared with individual gene markers in internal validation. Pending external validation in independent cohorts, it may provide molecular information that complements conventional histopathological evaluation in GC.

Key Words: Helicobacter pylori; Gastric cancer; Virulence factors; Transcriptional profiling; Gene expression; Stage-associated transcriptional signature; Disease staging

Core Tip: Standard diagnostic approaches for Helicobacter pylori generally establish whether the organism is present but do not resolve heterogeneity within gastric cancer. This study evaluates an expression-based discrimination model using a combined four-gene virulence transcript signature (cytotoxin-associated gene A, vacuolating cytotoxin A, blood group antigen-binding adhesin A, and induced by contact with epithelium gene A). The model identifies a stage-associated “adherence-to-invasion” transcriptional shift and, through a decision-tree classifier, discriminates early from advanced clinical stages with an area under the curve of 0.867, falling to 0.68 on nested cross-validation, a promising but modest signal warranting external validation.

INTRODUCTION

Gastric cancer (GC) continues to be among the leading causes of cancer-related deaths worldwide, with clinical outcomes strongly determined by stage at diagnosis and by the capacity to assess disease progression[1]. Helicobacter pylori(H. pylori) has a central role in gastric carcinogenesis; nevertheless, existing diagnostic strategies do not adequately capture the biological diversity that may underlie distinct clinical trajectories.

Efforts to improve risk assessment have largely centered on established virulence factors, especially cytotoxin-associated gene A (cagA) and vacuolating cytotoxin A (vacA), because of their strong associations with cancer development[2]. Such approaches remain limited, however, because they emphasize gene presence rather than functional activity and consequently yield restricted prognostic information. Accumulating evidence suggests that disease progression is shaped not simply by which genes are present, but also by the biological context in which virulence factors are expressed and interact[3].

Evidence from previous studies shows that H. pylori gene expression within the gastric tumor microenvironment can modify host pathways[4], immune responses[5], and microbial communities, each of which evolves dynamically across disease stages[6,7]. Most studies to date have compared the presence or absence of virulence genes among clinical groups such as gastritis, peptic ulcer disease, and cancer, rather than examining expression behavior after GC has become established. Consequently, an important stage-related knowledge gap remains: Integrated clinical models capable of translating bacterial activity into information relevant to treatment and prognosis, tailored to GC stage, are still lacking.

We hypothesize that stage-dependent variation in the expression of key H. pylori virulence genes reflects distinct functional disease states and can be modeled to estimate disease stage and clinical risk. We therefore developed a two-level molecular framework linking H. pylori virulence architecture and transcriptional activity with clinical disease stage. Population-based screening was integrated with detailed expression profiling of cagA, vacA, blood group antigen-binding adhesin A (babA), and induced by contact with epithelium gene A (iceA), followed by construction of a multivariate model designed to capture microbial expression patterns across stages and predict clinical outcomes. This approach shifts the analytical focus from pathogen identification toward transcriptional characterization, allowing transcriptional signatures to be considered as potential biomarkers of GC stage. The work extends a broader molecular characterization effort in this population; in the companion analysis, unsupervised clustering was used to jointly evaluate virulence and resistance characteristics across the entire cohort[8]. Here, the analysis is specifically directed toward the transcriptional trajectory associated with disease stage in the H. pylori-positive sub-cohort.

We term this pattern of declining babA expression alongside relatively increased cagA/vacA expression in advanced stages an “adherence-to-invasion” transcriptional shift - an expression pattern rather than a demonstrated biological transition. Accordingly, the present study treats H. pylori as a dynamic biomarker system whose behavior varies across disease stages rather than only as a static etiologic factor. Among patients with GC, the expression-based discrimination model may help identify individuals more likely to have advanced (stages III-IV) instead of early (stages I-II) disease at biopsy and may thereby inform management and referral decisions.

MATERIALS AND METHODS

Study design and stratification

This study employed a two-tiered design combining a case-control component (Tier one: 200 GC cases vs 200 controls) with a cross-sectional component (Tier two: Stage-associated expression analysis within the H. pylori-positive case sub-cohort), obtained from a previously established cohort. A subset of the same cohort had previously been evaluated using unsupervised clustering to investigate combined virulence and antimicrobial-resistance features[8]. For the present work, a two-tier analytical strategy was used to characterize the relationship between H. pylori transcriptional profiles and clinical stage. The primary cohort comprised 400 participants recruited at tertiary care centers in Pakistan: 200 individuals with histologically confirmed GC and 200 non-malignant controls.

For transcriptional analysis and predictive modeling, a sub-cohort of 159 H. pylori-positive GC cases was defined. Clinical stage was assigned according to the American Joint Committee on Cancer tumor-node-metastasis (TNM) classification system, specifically the 8th edition[9]. Subsequent analyses also grouped cases as early disease (stages I-II) or advanced disease (stages III-IV).

Tier one: Comparative virulence architecture (population level)

Polymerase chain reaction (PCR) was performed according to established protocols[10] to assess the prevalence of the H. pylori virulence genes cagA, vacA, babA, and iceA. Infection status was determined using a composite reference standard that integrated histological evaluation and molecular detection; a sample was considered positive if amplification was obtained for at least one target gene. Infection status was corroborated using species-conserved housekeeping genes glmM, ureA, hsp60 in addition to the virulence genes, thus reducing reliance on virulence-gene amplification alone.

The relationships between individual virulence markers and GC were quantified as odds ratios (ORs) with 95%CIs. Overall virulence burden was characterized by grouping samples by the number of detected virulence genes (X1-X4) and comparing these distributions across clinical stages. Supplementary material Sections 1 and 2 provide details of primer design, validation, and PCR conditions.

Tier two: Quantitative transcript profiling of H. pylori virulence genes (sub-cohort analysis)

Study population and cohort derivation: A two-tiered analytical design was employed. The primary study population comprised 200 cases of histologically confirmed gastric malignancy and 200 non-malignant controls, recruited from collaborating tertiary care centers in Pakistan. Cases were sourced from three tertiary care hospitals in Islamabad and Rawalpindi: Pakistan Institute of Medical Sciences, Holy Family Hospital, and the Pak Emirates Military Hospital in Rawalpindi[8]. Molecular analyses were performed at the Cancer Genetics and Epigenetics Laboratory, COMSATS University Islamabad (CUI). Specimens consisted of archived biopsies obtained retrospectively from an established tissue repository. The CUI Ethics Review Board approved the use of this repository in October 2017 (reference CIIT/Bio/ERB/17/15), and participant enrollment continued through 2021.

Primary study cohort (n = 200): The primary cohort included GC cases with histological confirmation of gastric malignancy and available clinicopathological information, including TNM stage. The control group comprised non-malignant individuals who underwent clinical evaluation with no histopathological evidence of gastric malignancy. Histological diagnoses were standardized according to the World Health Organization classification. All subtypes of adenocarcinoma were defined as adenocarcinoma, whereas non-adenocarcinoma gastric malignancies, such as neuroendocrine tumors, lymphomas, and gastrointestinal stromal tumors, were classified separately, in alignment with the companion cohort analysis[8]. Of the 200 cases in the primary cohort, 129 (64.5%) were adenocarcinoma and 71 (35.5%) were non-adenocarcinoma gastric malignancies. Among the 159 H. pylori-positive Tier Two subjects, 104 (65.4%) were adenocarcinoma and 55 (34.6%) were non-adenocarcinoma gastric malignancies.

Derivation of the H. pylori-positive analytical cohort (n = 159): Following initial molecular screening, all GC specimens underwent PCR-based assessment for H. pylori infection by detecting virulence-associated genes (cagA, vacA, babA, and iceA) and three species-conserved housekeeping genes (glmM, ureA, and hsp60). Among the 200 GC cases, 159 (79.5%) had at least one detectable H. pylori virulence marker and/or housekeeping gene and were therefore eligible for Tier Two transcriptomic analysis. The other 41 cases were omitted from gene expression analyses because no amplification of any virulence or housekeeping target was detected. Thus, quantitative transcriptional profiling, multivariate modeling, and stage-stratified analyses were all based on the 159 GC specimens positive for H. pylori. Supplementary Figure 1 summarizes this participant flow in diagram form.

Inclusion and exclusion criteria: Eligibility for analysis required fulfillment of four criteria: (1) Histologically confirmed gastric malignancy (adenocarcinoma or other GC histology); (2) Gastric biopsy tissue available and suitable for DNA and RNA extraction; (3) H. pylori genomic DNA detectable by PCR; and (4) Nucleic acid quality adequate for downstream molecular analyses.

Specimens were excluded because of inadequate tissue quantity, poor nucleic acid quality, severe degradation or chemical cross-linking, previous neoadjuvant chemotherapy or radiotherapy, prior total gastrectomy, or serious gastric comorbidities considered likely to compromise the reliability of tissue and molecular data. Complete inclusion and exclusion details are available in Supplementary material Section 3.

Sample size justification: Archived gastric biopsy specimens from an established GC cohort formed the basis of this analysis. Every archived specimen satisfying the prespecified molecular and clinical eligibility criteria was included. After H. pylori screening, clinicopathological stratification, and nucleic acid quality control, 159 H. pylori-positive GC specimens remained available for downstream transcriptomic analysis.

Adequacy of the analytical cohort was examined using a power-based sample size calculation informed by assumptions from earlier expression studies and pilot observations. Assuming a standard deviation of 0.5, an expected mean difference of 0.3 in relative mRNA expression, a two-sided significance level of 0.05, and 80% statistical power, the calculation indicated a minimum requirement of 44 participants per group for expression-based comparisons[11]. Stage III (n = 18) and IV (n = 37) in the H. pylori-positive sub-cohort fell short of this per-group target (see limitation paragraph)].

RNA extraction and complementary DNA synthesis: Total RNA was isolated from gastric biopsy tissue by a TRIzol™-based dual-phase extraction procedure that also permitted simultaneous recovery of genomic DNA. RNA purity and integrity were assessed by NanoDrop™ spectrophotometry and 1% agarose gel electrophoresis. Further analysis was restricted to samples that satisfied the prespecified quality criteria, including the A260/280 ratio[12,13]. Complementary DNA was generated with the SuperScript™ III First-Strand Synthesis System (Invitrogen) and random hexamer primers under standardized reaction conditions (Supplementary material Section 2).

Quantitative real-time PCR: SYBR™ Green-based real-time PCR on the ABI PRISM 7000/7500 platform was used to quantify expression of cagA, vacA, babA, and iceA. Reactions were run in triplicate, and the mean Ct value from each set of replicates was used for analysis. To limit normalization bias arising from differences in bacterial load and tissue quality, relative expression was normalized to the geometric mean of three housekeeping genes consistently expressed in H. pylori (ureA, glmM, and hsp60). Amplifiable cDNA was obtained from all 159 eligible specimens, which were consequently retained for subsequent analyses.

Expression fold changes relative to the control cohort were calculated with the 2^-ΔΔCt method. Supplementary material Sections 1 and 2 contain the complete primer sequences, in silico validation procedures, and reaction conditions. Primers were designed and bioinformatically validated using NCBI Primer-BLAST, PCR Primer Stats, and the Sequence Manipulation Suite (Supplementary material Section 1.4). Assay specificity was assessed by melt-curve analysis and triplicate reactions for every sample. Amplification efficiency was verified during assay optimization per standard quantitative PCR (qPCR) practice. No-template and no-reverse-transcription controls were included for every target and housekeeping gene assay to exclude reagent contamination and genomic DNA carryover, respectively.

Computational modeling and stage-associated expression analysis

A multi-parametric analytical framework was implemented in the R environment (v4.5.2) to translate complex molecular expression measurements into predictive models for stage classification. Regression analyses used the MASS package; recursive partitioning used rpart and rpart.plot; principal component analysis (PCA) used stats; expression-based discrimination models were evaluated with pROC; and robust model validation was conducted with caret.

Ordinal logistic regression: Ordinal logistic regression (OLR) models examined associations between virulence gene expression and disease stage, with TNM stage specified as an ordered categorical outcome. Model coefficients, ORs, and 95%CIs were estimated, and model fit was assessed with likelihood ratio statistics and Akaike information criterion (AIC) values. The proportional-odds assumption was tested for each predictor using the Brant test. Where violated, the reported OR is interpreted as an average effect across stage cutpoints rather than a uniform per-cutpoint effect (see limitations paragraph).

Classification and regression tree: Cases were classified as early (stages I-II) or advanced (stages III-IV) disease using recursive partitioning with the classification and regression tree (CART) algorithm. Gini impurity served as the principal splitting criterion. Tree complexity was optimized using cost-complexity pruning, with the complexity parameter selected based on minimum cross-validation error. Tuning and pruning were nested within an outer 5-fold cross-validation, so no held-out fold influenced its own tree's tuning. Performance metrics were pooled across the five held-out folds.

Principal component analysis: Principal component analysis was conducted on mean-centered, unit-variance-scaled expression values to evaluate the multivariate structure of the gene expression data and to visualize sample clustering across clinical stages.

Model validation and robustness testing: A multimetric validation strategy was used to examine the clinical reliability and generalizability of the predictive framework. Discriminatory performance was quantified by receiver operating characteristic (ROC) analysis, using area under the curve (AUC) as the primary metric. Separate ROC analyses evaluated the multivariate ordinal regression framework and the binary classifier differentiating early-stage (stages I-II) from advanced-stage (stages III-IV) cases. For the ordinal regression framework, the fitted linear predictor combining all four gene-expression values served as a continuous risk score for this early/advanced-stage discrimination task. Sensitivity and specificity were additionally examined to determine clinically relevant operating thresholds.

Model stability and potential overfitting were evaluated through 5-fold cross-validation implemented in the caret framework. Robustness of the OLR model was additionally examined with a Cauchit link specification. For the CART classifier, tuning and evaluation were nested within an outer 5-fold cross-validation, so that the reported AUC, sensitivity, and specificity reflect performance on held-out folds rather than the training data. The ordinal regression-based combined signature was evaluated analogously: The four-gene model was refit on each training fold, and the resulting linear predictor was applied to the corresponding held-out fold, with AUC and calibration pooled across all five held-out folds. Variance inflation factors (VIFs) were calculated to assess multicollinearity among predictors; values below 5 were interpreted as evidence of acceptable predictor independence.

All statistical tests were two-sided unless otherwise specified, with P < 0.05 considered statistically significant. All analytic variables in both Tiers had complete data. Given the exploratory four-gene panel evaluated across multiple stages, models, and pairwise correlations, P values were not corrected for multiplicity. Accordingly, these findings should be interpreted as hypothesis-generating and require independent replication.

RESULTS

Development of a stage-associated transcriptional framework for GC

We developed and evaluated a framework for identifying H. pylori virulence transcriptional signatures associated with GC stage. H. pylori virulence markers and healthy controls provided the molecular reference framework. Tier one consisted of population-level molecular screening designed to separate malignant cases from healthy controls; Tier two used quantitative transcript modeling in H. pylori-positive patients to examine transcriptional patterns in relation to clinical stage.

Tier one: Clinical characteristics and population-level molecular screening

Tier one established the molecular basis for initial GC screening. Across the 400 participants, tumor cases and healthy controls were compared for prevalence and disease associations of the four-gene virulence panel (cagA, vacA, babA, iceA). The analysis identified the factors with the strongest relationships to GC and provided the foundation for the subsequent Tier 2 assessment of stage-associated transcriptional variation.

Baseline demographics and pathological characteristics: The primary set of 200 GC cases was analyzed to determine whether host demographic or pathological characteristics varied with disease progression. Mean age was 52.9 ± 16.4 years, and male participants accounted for 56.0% of the group. Neither age nor sex was significantly associated with TNM stage (P = 0.532 and P = 0.727, respectively). Likewise, the lifestyle and anatomical characteristics reported, including smoking status (P = 0.770) and tumor anatomical location (P = 0.173), showed no significant relationship with disease progression. Distal (28.5%) and proximal (27.5%) tumors were frequent, but primary tumor site remained independent of clinical stage (Figure 1-C2).

Figure 1 Clinicopathological and molecular profile of gastric cancer (n = 200). A: Clinical progression (A1 and A2) demonstrates a stage-dependent movement toward poorer tumor differentiation (P = 0.001), accompanied by a pronounced reduction in cure rates from 55.9% in stage I to 4.4% in stage IV (P < 0.001); B: Virulence-associated risk (B1-B2) compares the distributions of Helicobacter pylori (H. pylori) genotypes in cases and controls. The strongest associations with malignancy are observed for the toxigenic factors, particularly vacuolating cytotoxin A [odds ratio (OR) = 22.9] and cytotoxin-associated gene A (OR = 18.7) (P < 0.001); C: Genomic complexity and anatomy. Panel C1 depicts stage-related differences in virulence-gene accumulation (P = 0.045): The quadruple genotype (X4) reaches its highest frequency in stage II, whereas the triple genotype (X3) occurs more frequently in advanced disease. Panel C2 presents the distribution of tumor anatomical location across four categories. Data shown in B1, B2, and C1 were obtained from the molecularly confirmed H. pylori-positive sub-cohort (n = 159) for identification of pathogen-specific risk factors. TNM: Tumor-node-metastasis; OR: Odds ratio; CI: Confidence interval.

Histopathological characteristics, by contrast, showed stage-related distributions. Histological diagnosis (adenocarcinoma vs other gastric malignancy histology) was significantly associated with stage (P = 0.034; Table 1), primarily due to a higher proportion of adenocarcinoma among Stage III cases (90.9%) compared to other stages (58.5%-63.2%). Tumor differentiation grade was associated with cancer stage (P = 0.001). Figure 1-A1 shows that stage I included substantial proportions of atypical (47.1%) and well-differentiated (20.6%) tumors. As disease advanced, the distribution shifted progressively toward poorer differentiation, with poorly differentiated tumors reaching a maximum of 45.5% in stage III. The decline in histopathological differentiation was accompanied by poorer survival outcomes: Cure rates decreased markedly from 55.9% in stage I to 4.4% in stage IV (P < 0.001; Figure 1-A2).

Table 1 Baseline demographics and significant clinical characteristics (n = 200), n (%).

Variable	Overall (n = 200)	Stage I (n = 68)	Stage II (n = 65)	Stage III (n = 22)	Stage IV (n = 45)	P value
Age (years)	52.9 ± 16.4	51.2 ± 15.4	54.3 ± 19.3	56.1 ± 12.4	52.0 ± 15.0	0.532
Gender						0.727
Female	88 (44.0)	27 (39.7)	32 (49.2)	10 (45.5)	19 (42.2)
Male	112 (56.0)	41 (60.3)	33 (50.8)	12 (54.5)	26 (57.8)
Histological Dx						0.034a
Adenocarcinoma	129 (64.5)	43 (63.2)	38 (58.5)	20 (90.9)	28 (62.2)
Others	71 (35.5)	25 (36.8)	27 (41.5)	2 (9.1)	17 (37.8)
Histological grade						0.001b
Well	47 (23.5)	14 (20.6)	18 (27.7)	6 (27.3)	9 (20.0)
Moderate	48 (24.0)	13 (19.1)	14 (21.5)	5 (22.7)	16 (35.6)
Poorly	47 (23.5)	9 (13.2)	14 (21.5)	10 (45.5)	14 (31.1)
Atypical	58 (29.0)	32 (47.1)	19 (29.2)	1 (4.5)	6 (13.3)
H. pylori (histology)						0.042a
Positive	130 (65.0)	36 (52.9)	48 (73.8)	17 (77.3)	29 (64.4)
Negative	70 (35.0)	32 (47.1)	17 (26.2)	5 (22.7)	16 (35.6)
Survival status						< 0.001b
Deceased	65 (32.5)	15 (22.1)	17 (26.2)	10 (45.5)	23 (51.1)
Under treatment	82 (41.0)	15 (22.1)	36 (55.4)	11 (50.0)	20 (44.4)
Cured	53 (26.5)	38 (55.9)	12 (18.5)	1 (4.5)	2 (4.4)

^aP < 0.05.

^bP < 0.01. One-way ANOVA was used to calculate P values for the continuous variable (age), whereas the χ² test or Fisher’s exact test, as appropriate, was used for categorical variables. H. pylori: Helicobacter pylori.

Comparative virulence prevalence and risk assessment: Molecular profiling characterized the prevalence of virulence factors throughout the study population. Histology detected H. pylori in 65.0% of cancer cases, and positivity was more frequent at advanced stages (P = 0.042). Molecular testing subsequently detected H. pylori in 79.5% (n = 159) of the cancer group. Each of the four H. pylori virulence genes evaluated (cagA, vacA, babA, and iceA) was enriched among GC cases.

Cancer risk associated with the bacterial virulence factors was then evaluated within the molecularly confirmed H. pylori-positive sub-cohort (n = 159). The resulting estimates showed disease associations that varied by gene.

Within this sub-cohort, virulence-associated cancer risk followed a clear statistical hierarchy (Figure 1-B1 and B2; Table 2). The strongest associations with GC were observed for the effector toxins vacA (OR = 22.9; 95%CI: 11.6-45.2) and cagA (OR = 18.7; 95%CI: 10.9-32.3) (P < 0.001). A significant risk association was also observed for the adhesion gene babA (OR = 14.1; 95%CI: 8.6-23.2), consistent with an important role for epithelial attachment in gastric carcinogenesis. iceA expression showed a smaller association (OR = 3.01; 95%CI: 1.84-4.94), although statistical significance remained high (P < 0.001).

Table 2 Prevalence and association of Helicobacter pylori virulence markers, n (%).

Marker	Gastric cancer positive (n = 159)	Controls (n = 200)	Odds ratio	95%CI	P value
cagA	131 (82.4)	40 (20.0)	18.7	10.9-32.3	< 0.001b
vacA	148 (93.1)	74 (37.0)	22.9	11.6-45.2	< 0.001b
babA	117 (73.6)	33 (16.5)	14.1	8.6-23.2	< 0.001b
iceA	58 (36.5)	32 (16.0)	3.01	1.84-4.94	< 0.001b

^bP < 0.01.

cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A.

Distribution of virulence factor accumulation by clinical stage: Virulence factor accumulation was also evaluated according to clinical stage. Figure 1-C1 and Table 3 show that multigene genotype frequencies varied across stages, with the triple (X3) and quadruple (X4) genotypes reaching their highest frequencies in later disease stages. The X4 genotype, containing all four genes, was especially frequent in stage II, where it reached 22.2% (P = 0.045). Overall, patterns of virulence-gene co-occurrence differed by stage: vacA and cagA remained highly prevalent, while babA and iceA displayed more variable stage-specific distributions.

Table 3 Distribution of virulence factor accumulation by clinical stage, n (%).

Virulence factor accumulation	Stage I (n = 50)	Stage II (n = 54)	Stage III (n = 18)	Stage IV (n = 37)	P value
1 gene	5 (10.0)	2 (3.7)	0 (0.0)	5 (13.5)
2 genes	5 (10.0)	4 (7.4)	5 (27.8)	8 (21.6)
3 genes	34 (68.0)	36 (66.7)	13 (72.2)	19 (51.4)
4 genes	6 (12.0)	12 (22.2)	0 (0.0)	5 (13.5)	0.045a

^aP < 0.05.

The “adherence-to-invasion” transcriptional pattern

Progression of GC was accompanied by a marked change in the H. pylori transcriptional profile. Specifically, the relative expression of adherence-associated genes decreased as disease progressed, whereas that of cytotoxin-related genes increased.

Stage-specific transcriptional remodeling of the virulence program: Joint evaluation of fold-induction kinetics and absolute transcript abundance demonstrated stage-dependent differences in virulence gene expression among the 159 infected cases. Fold-induction profiles indicated a transition from an adherence-oriented pattern in early disease to a toxin-oriented phenotype with progression (Table 4; Figure 2A). The distributions of absolute mRNA levels independently supported the same stage-specific pattern (Table 5; Figure 2B).

Figure 2 Transcriptional transition from bacterial adherence to host-cell invasion across clinical stages of gastric cancer. A: Fold induction of major virulence transcripts relative to baseline expression in asymptomatic controls is shown, with the dashed horizontal line indicating y = 1. A distinct “crossover” is evident at stage II: Expression of the adhesin gene blood group antigen-binding adhesin A (babA) decreases, while cytotoxin-associated gene A (cagA) and vacuolating cytotoxin A (vacA), the toxigenic factors, become progressively upregulated; B: Relative mRNA expression of cagA, vacA, babA, and induced by contact with epithelium gene A according to stage. Values are expressed as 2^-ΔΔCt; boxes indicate the interquartile range (IQR), horizontal lines indicate medians, and whiskers extend to 1.5 × IQR; C: The Spearman’s rank correlation matrix summarizes transcriptional relationships among virulence factors. Correlation coefficients (r) are represented by the color scale, including the significant positive relationship between cagA and vacA (r = 0.31); D: Principal component analysis of molecular profiles demonstrates progressive directional separation with 90% confidence ellipses, consistent with a global microbial-transcriptome shift from early-stage (stage I) to advanced (stage IV) gastric cancer. PCA: Principal component analysis; cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A.

Table 4 Fold induction of Helicobacter pylori virulence genes relative to healthy controls (n = 159).

Clinical stage	cagA	vacA	babA	iceA
Stage I	3.11	3.15	5.68	0.83
Stage II	4.81	2.90	1.80	2.35
Stage III	4.53	4.06	0.67	0.68
Stage IV	4.26	2.41	1.59	0.52

cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A.

Table 5 Relative mRNA expression of Helicobacter pylori virulence factors by tumor-node-metastasis stage (n = 159).

Virulence gene	Stage I median (IQR)	Stage II median (IQR)	Stage III median (IQR)	Stage IV median (IQR)	P value
cagA	1.93 (1.04-3.36)	3.20 (1.51-5.98)	3.80 (1.57-6.49)	3.11 (1.09-6.20)	< 0.001b
vacA	2.06 (1.36-3.48)	2.42 (1.41-3.92)	3.13 (1.85-3.91)	1.79 (1.12-3.20)	< 0.001b
babA	4.11 (1.62-5.41)	1.39 (1.03-2.11)	0.57 (0.30-1.11)	1.19 (0.54-1.72)	< 0.001b
iceA	0.63 (0.39-0.94)	0.86 (0.46-3.17)	0.54 (0.40-0.91)	0.65 (0.30-1.41)	0.027a

^aP < 0.05.

^bP < 0.01; the Kruskal-Wallis test was used for comparisons across the four stages. cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A; IQR: Interquartile range.

The oncogenic effector cagA increased significantly from a median of 1.93 [interquartile range (IQR) 1.04-3.36] in stage I to its maximum median value of 3.80 (IQR 1.57-6.49) in stage III (P < 0.001). vacA likewise reached its highest expression in stage III, with a median of 3.13 (IQR 1.85-3.91; P < 0.001). Together, these trajectories are consistent with progressive strengthening of the toxigenic axis as clinical stage advances.

The opposite pattern was observed for the adhesin gene babA. Expression was highest during early disease, with a median of 4.11 (IQR 1.62-5.41) in stage I, but fell sharply to 0.57 (IQR 0.30-1.11) by stage III (P < 0.001), corresponding to an approximately sevenfold reduction in transcriptional output. iceA transcript levels remained relatively low throughout the stages; a modest but significant rise occurred in stage II (median 0.86; IQR 0.46-3.17; P = 0.027), followed by lower expression at later stages. Taken together, the data describe a stage-associated reconfiguration of virulence expression in which adherence-related transcription is more prominent in early disease and toxigenic activity becomes more prominent in advanced disease (Figure 2A and B).

A statistically significant moderate positive correlation between cagA and vacA was observed in Table 6 and Figure 2C (r = 0.31, P < 0.001), consistent with coordinated transcription within the principal toxigenic axis. Relationships involving the secondary factors were weaker. babA had a modest positive correlation with vacA (r = 0.22), whereas iceA was mildly inversely correlated with cagA (r = -0.18).

Table 6 Spearman correlation matrix of virulence factor expression.

Gene	cagA	vacA	babA	iceA
cagA	1.00	0.31	0.08	-0.18
vacA	0.31	1.00	0.22	-0.06
babA	0.08	0.22	1.00	-0.03
iceA	-0.18	-0.06	-0.03	1.00

cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A.

PCA was used to characterize the broader organization of the transcriptional network (Figure 2D). Stage groups showed differing distributions in multivariate expression space, although the 90% confidence ellipses for adjacent stages overlapped substantially, indicating partial rather than complete stage-associated separation of overall virulence expression profiles. Considered together, the correlation map (Table 6; Figure 2C) and PCA results (Figure 2D) indicate a coordinated H. pylori virulence-expression network dominated by the cagA-vacA module and further shaped by stage-dependent changes in adherence-related factors.

Modeling of associations with disease stage

The association between H. pylori transcriptional activity and clinical stage was evaluated by fitting an OLR model in the H. pylori-positive group (n = 159). In the multivariate model, babA and iceA were the most significant independent transcriptional factors related to clinical stage.

Among the modeled transcripts, reduced babA expression had the strongest relationship with advanced disease (OR = 0.536, 95%CI: 0.428-0.654, P < 0.001). iceA also retained an independent significant association (OR = 0.802, P = 0.010). By comparison, the canonical toxigenic markers cagA (P = 0.123) and vacA (P = 0.924) did not independently reach statistical significance in the multivariate analysis (Figure 3A, Table 7).

Figure 3 Validation and discriminative utility of the multigene transcriptional signature (n = 159). A: Odds ratios (OR) with 95%CIs are displayed. Blood group antigen-binding adhesin A (babA) (OR = 0.536) and induced by contact with epithelium gene A (OR = 0.802) are the only independent discriminators of advanced-stage disease that reach statistical significance (purple markers, P < 0.05); B: Comparison of Akaike information criterion scores indicates superior fit for the full panel (377) relative to the univariate (423.3) and dual-toxin (424.3) models; C: Receiver operating characteristic analysis shows that the combined signature (thick black line) reaches an area under the curve of 0.769 (apparent; falling to 0.747 on pooled 5-fold cross-validation) and performs better than the individual markers; D: Boxplots show a significant relationship between pathological stage and transcriptional scores (Spearman r = 0.494, P < 0.001), with babA downregulation providing the strongest association. ROC: Receiver operating characteristic; cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A.

Table 7 Multivariate ordinal logistic regression parameter estimates (n = 159).

Predictor/threshold	(β)	Std. Error	z-statistic	P value	OR	95%CI
babA	-0.624	0.107	-5.814	< 0.001b	0.536	0.428-0.654
iceA	-0.221	0.085	-2.585	0.010a	0.802	0.675-0.945
cagA	0.073	0.047	1.544	0.123	1.076	0.981-1.182
vacA	0.009	0.091	0.095	0.924	1.009	0.843-1.205
Stage I|II	-2.368	0.441	-5.376	< 0.001b	-	-
Stage II|III	-0.472	0.393	-1.201	0.230	-	-
Stage III|IV	0.176	0.389	0.452	0.651	-	-

^aP < 0.05.

^bP < 0.01.

cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A; OR: Odds ratio.

Model selection and statistical optimization: Three candidate models were compared to identify the most robust approach for discriminating disease stage: The univariate model (cagA), the dual-toxin model (cagA + vacA), and the complete multivariate model (babA + iceA + cagA + vacA) (Table 8).

Table 8 Model selection and comparison statistics.

Model specification	Log-likelihood	AIC (lower is better)	Model χ2	P value
Univariate (cagA)	-207.67	423.33	3.30	0.069
Toxigenic axis (cagA + vacA)	-207.17	424.34	4.29	0.117
Full multivariate panel	-181.50	377.01	55.63	< 0.001b

^bP < 0.01.

AIC: Akaike information criterion; cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A.

Relative to the null baseline, the univariate model offered little explanatory improvement (χ² = 3.30, P = 0.069). Adding a second toxin marker yielded only a small improvement in model fit, whereas the complete multivariate panel performed best. Incorporation of all four virulence genes significantly improved fit (χ² = 55.63, P < 0.001) and lowered the AIC to 377.01 (Figure 3B). Consistent with these results, the likelihood ratio test favored the four-gene panel for stage-discriminative transcriptional modeling (P < 0.001).

Multivariate parameter significance: OLR threshold estimates indicated that the stage I-to-stage II transition (intercept = -2.368, P < 0.001) accounted for the largest transcriptional change captured by the model. Separation was less distinct for the stage II-to-stage III and stage III-to-stage IV transitions, implying that the strongest transcriptional differentiation occurs relatively early in disease stratification.

The performance of the integrated signature was further characterized using standard classification measures, including sensitivity and specificity across disease stages. To reduce potential confounding, we refitted the four-gene model and adjusted for age, sex, smoking status, tumor location, and differentiation grade. After adjustment, babA (adjusted OR = 0.524, 95%CI: 0.420-0.653, P < 0.001) and iceA (adjusted OR = 0.796, 95%CI: 0.663-0.955, P = 0.014) remained significant independent predictors of stage. Sex and smoking status were also independently linked to stage in this model.

To assess the added value of the transcriptional panel compared to routine clinical variables, we compared the model with only age, sex, smoking status, tumor location, and differentiation grade to a combined clinical-plus-gene model using pooled 5-fold cross-validated AUC for early-versus-advanced-stage discrimination. The clinical-only model performed at chance level (AUC = 0.49, 95%CI: 0.40-0.58). Adding the four-gene panel increased discrimination to AUC = 0.71 (95%CI: 0.62-0.79; bootstrap 95%CI for the difference: 0.12-0.32). The four-gene signature alone performed comparably or better (AUC = 0.75) than the combined clinical-plus-gene model (Figure 4A).

Figure 4 Discrimination, incremental value, and calibration of the combined four-gene transcriptional signature (n = 159). A: Discrimination across model specifications was evaluated using receiver operating characteristic area under the curve (AUC) for distinguishing early (stages I-II) from advanced (stages III-IV) cases. The four-gene panel alone achieved AUC = 0.747 (cross-validated), vs AUC = 0.49 for the clinical-variables-only model. The dotted reference line represents chance-level discrimination (AUC = 0.5); B: Calibration assessed by comparing the observed proportion of advanced-stage (III-IV) cases to the mean model-predicted probability across five equal-sized bins of pooled out-of-fold predicted probabilities. The dashed diagonal line indicates perfect calibration. Calibration was adequate (Hosmer-Lemeshow χ² = 5.53, df = 3, P > 0.05). AUC: Area under the curve; ROC: Receiver operating characteristic; CV: Cross-validation.

Model validation and stage-discrimination performance: The integrated molecular signature achieved an AUC of 0.769 and outperformed every individual gene predictor (Table 9; Figure 3C). This apparent AUC fell modestly to a pooled 5-fold cross-validated AUC of 0.747 (95%CI: 0.664-0.823), indicating good generalizability; calibration of the cross-validated predicted probabilities was adequate (Hosmer-Lemeshow χ² = 5.53, df = 3, P > 0.05), with no evidence of systematic over- or under-prediction across risk strata (Figure 4). VIF values clustered near 1 (mean ≈ 1.05), demonstrating minimal multicollinearity among the four transcripts and supporting the interpretation that each transcript contributed largely independent information to the model.

Table 9 Stage-associated transcriptional signature validation and correlation analysis (n = 159).

Parameter/variable	Spearman (ρ)	P value	VIF	ROC-AUC (95%CI)	Directional trend/validation
babA	-0.467	< 0.001b	1.05	0.702 (0.632-0.795)	Downregulation
iceA	-0.133	0.155	1.00	0.584 (0.540-0.717)	Variable
cagA	0.122	0.126	1.05	0.612 (0.453-0.645)	Upregulation
vacA	-0.010	0.210	1.10	0.528 (0.457-0.647)	Neutral
Combined model (apparent, full-cohort)	ρ = 0.494	< 0.001b	1.02	0.769	Linear fit
Combined model (Nested 5-fold CV)	-	-	-	0.747 (0.664-0.823)	Calibration: Hosmer-Lemeshow χ2 = 5.53, df = 3, P > 0.05
5-fold cross-validation (Cauchit)	-	-	-	-	Accuracy: 0.535; Kappa: 0.325
CART - decision tree (apparent, full-cohort)	-	-	-	0.867	Sensitivity: 74.5%; specificity: 81.7%
CART - decision tree (Nested 5-fold CV)	-	-	-	0.68 (95%CI: 0.59-0.77)	Sensitivity: 54.5%; specificity: 75.0%

^bP < 0.01.

cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A; VIF: Variance inflation factor; ROC: Receiver operating characteristic; AUC: Area under the curve; CART: Classification and regression tree; CV: Cross-validation.

A significant Spearman correlation between predicted and observed stages confirmed concordance between the two measures (r_s = 0.494, P < 0.001) (Figure 3D). Cross-validation also supported internal consistency. With a Cauchit link function, five-fold cross-validation produced a mean accuracy of 0.535 and a Kappa coefficient of 0.325. Although these values are modest, they indicate performance above chance. They should nevertheless be interpreted as resampled internal performance within this cohort, rather than as evidence of external validity or immediate clinical applicability.

Because Stage III contained fewer cases, a sensitivity analysis was performed after merging Stage III and Stage IV into a single category. Refitting the same Cauchit-link model to the combined outcome yielded performance comparable with or superior to that of the four-stage analysis (5-fold cross-validation accuracy = 0.604, Kappa = 0.403). Both babA (P < 0.001) and iceA (P = 0.004) again remained the only statistically significant predictors, and their coefficients changed minimally relative to the four-stage model.

The four-gene transcriptional signature was next translated into a practical framework for differentiating early from advanced disease. Recursive partitioning was used to construct a decision-tree model that classified patients into early (stages I-II) and advanced (stages III-IV) risk groups, while identifying a transcriptional threshold associated with advanced-stage disease. Cost-complexity pruning, combined with cross-validation, was applied to control tree complexity, promote generalizability, and reduce overfitting.

In Figure 5A, the initial and most informative split was defined by babA expression at a threshold ≥ 0.72, making babA the principal sentinel marker for separation of early-stage and advanced-stage disease. Additional splits involving cagA, iceA, and vacA refined classification within the higher-risk groups. ROC analysis of the resulting stage-associated transcriptional signature yielded an AUC of 0.867. A nested 5-fold cross-validation, with tuning and pruning confined to each training fold, gave a more conservative AUC of 0.68 (95%CI: 0.59-0.77; sensitivity 54.5%, specificity 75.0%). We report the nested estimate as the classifier’s validated performance (Figure 5B). The apparent-versus-cross-validated comparison across all three metrics is summarized in Figure 5C.

Figure 5 Predictive model, performance validation, and apparent-versus-cross-validated comparison for stage groups. A: Recursive partitioning (decision tree) classifier separating patients into early (stages I-II) and advanced (stages III-IV) groups. Blood group antigen-binding adhesin A expression at a threshold ≥ 0.72 defines the first branch, after which cytotoxin-associated gene A, induced by contact with epithelium gene A, and vacuolating cytotoxin A provide additional downstream splits. Each terminal node reports the predicted class (early in blue, advanced in orange), the probability of class membership, and the proportion of the cohort; B: Receiver operating characteristic analysis evaluating stage discrimination by the integrated four-gene signature. The tree fit on the full cohort showed an apparent area under the curve (AUC) of 0.867; nested 5-fold cross-validation gave a more conservative AUC of 0.68. The diagonal line represents chance-level performance; C: Apparent vs nested cross-validated performance of the classification and regression tree classifier across three metrics. The consistent decrease from apparent to cross-validated estimates across all three metrics (apparent AUC = 0.867 vs nested CV AUC = 0.68) shows that the apparent fit is optimistic, as expected for an in-sample tree evaluation. cagA: Cytotoxin-associated gene A; vacA: Vacuolating cytotoxin A; babA: Blood group antigen-binding adhesin A; iceA: Induced by contact with epithelium gene A; AUC: Area under the curve; CV: Cross-validation.

DISCUSSION

The present study establishes a molecular framework that connects H. pylori virulence dynamics with progression of GC. Integration of population-level screening and transcriptional profiling identified systematic stage-associated differences in H. pylori virulence expression during tumor progression. A transcriptional threshold associated with stage was characterized by a shift from adherence-oriented expression toward toxin-oriented expression. These states may serve as functional molecular indicators of disease severity and may provide a quantitative adjunct to conventional histopathology.

H. pylori virulence as a multidimensional correlate of gastric carcinogenesis

Tier one demonstrated that the established virulence factors cagA and vacA were significantly more prevalent among GC cases than among healthy controls, further supporting the recognized importance of the cag pathogenicity island and vacA as major pathogenic elements in oncogenic processes[2,14]. These observations extend our previous characterization of virulence-associated molecular patterns in this population, in which oncogenic H. pylori genotypes were highly prevalent among patients with GC[8]. In the companion analysis, unsupervised clustering of combined virulence and antimicrobial-resistance features was used to define case-level molecular profiles across the full cohort. The current study addresses a different question by examining the stage-associated transcriptional trajectory within the H. pylori-positive sub-cohort and determining how effectively that trajectory discriminates TNM stages.

These results indicate that infection status by itself does not adequately predict clinical progression. Instead, distinct pathogenic states appear to be defined by the combined accumulation and coordinated expression of virulence determinants. This interpretation is consistent with the emerging view that multi-marker virulence profiles provide more informative prediction of adverse clinical outcomes than individual genetic markers[15].

Evidence of a stage-associated signal was already visible at the level of pathogen detection. Histology identified H. pylori in 65.0% of cases overall and showed a disproportionately higher detection frequency in advanced disease (P = 0.042). Molecular detection based on virulence-associated genes, however, identified infection in a substantially greater proportion of the cohort (79.5%). The difference may be explained partly by limited histological sensitivity, particularly when bacterial loads are low[16]. More broadly, these findings emphasize that assessment based solely on H. pylori presence may underestimate infection and fail to capture stage-related information revealed by molecular characterization.

Recent studies provide additional context for this multidimensional risk by showing that aggressive H. pylori genotypes can reshape the gastric microbial environment and facilitate cancer development[6,17]. Such interactions support a shift away from simple “positive/negative” testing toward a more detailed characterization of virulence activity states when assessing cancer risk across disease stages.

The composition of the reference group requires particular consideration. Table 2 shows that per-gene carriage among controls ranged from 16.0% to 37.0%, indicating that at least 37.0% of controls carried one or more of the four virulence genes assessed. As with case classification, designation of controls as “H. pylori-negative” was based solely on the seven-gene composite panel (cagA, vacA, babA, iceA, glmM, ureA, and hsp60), incorporating three species-conserved housekeeping genes alongside the four polymorphic virulence genes. This reduces, though does not entirely eliminate, the chance of misclassifying a truly infected individual as negative, as some controls classified as negative may still harbor H. pylori genes below this panel’s PCR detection threshold. Actual background H. pylori carriage in this population is likely higher still, as asymptomatic prevalence in Pakistan has been reported at 45%-74%, depending on region and detection method[18,19]. If undetected carriers were present among the controls, the odds ratios reported in Table 2 would represent risk relative to a reference group that was not completely H. pylori-free. Such misclassification would be expected to bias the associations toward the null rather than exaggerate them.

Transcriptional remodeling: The “adherence-to-invasion” transcriptional pattern

Tier two transcriptome analyses identified a distinct stage-dependent reorganization of virulence expression. As the disease progressed, expression of adherence factors such as babA decreased, while expression of cytotoxin-associated genes such as cagA and vacA increased. The pattern provides molecular evidence that disease severity is associated with specific transcriptional programs rather than with the mere carriage of individual genes. Correlation results further indicate that cagA and vacA operate as coordinated elements of the toxigenic system, whereas adherence regulators such as babA behave with greater relative independence. These contrasting patterns may reflect changing biological requirements within the bacterial niche as malignant transformation progresses.

Previous work has shown that BabA expression is more frequent in strains isolated from patients with peptic ulcers and GC than in strains associated with asymptomatic colonization[20]. Most of that literature compares separate clinical categories, including infected vs uninfected or symptomatic vs asymptomatic groups, rather than following variation within patients who already have established cancer. The present results address this different dimension: Among confirmed GC cases, babA expression decreases specifically as TNM stage advances. Existing BabA studies have only rarely examined this within-disease trajectory directly.

This divergence is consistent with the hypothesis that H. pylori transcription shifts from a colonization-associated profile toward an invasion-associated pattern as disease advances[21]. Bacterial invasion itself was not directly measured, so this interpretation remains speculative. The repression of babA from Stage I through III may reflect a molecular “detachment” pattern, possibly related to changing host microenvironments or immune pressures[20,22]. However, babA expression rose again at Stage IV (Table 5) - a non-monotonic pattern not explained by a simple linear adherence-to-invasion trajectory. Whether this reflects broader multi-omics divergence in the gastric tumor microbiome during carcinogenesis[23,24] was not tested directly in this cohort.

Mechanistic implications and biological interpretation

Chronic inflammation and bacterial signaling may plausibly contribute to this transcriptional pattern, though neither was measured directly here. One speculative interpretation is that, in a malignant environment, the combination of high cytotoxin expression and low adhesion could reflect an evolutionary adaptation[25]. Reduced epithelial anchoring could alter host-pathogen interactions, while persistent cytotoxin expression could maintain biological activity in the tumor microenvironment[2]. Recent evidence further suggests that microbial reprogramming may influence ferroptosis, immune checkpoint regulation, and formation of advanced tumors[26]. Whether bacterial virulence gene expression influences tumor cell metabolism or immune surveillance in this cohort, vs simply correlating with disease stage, was not tested here and remains speculative.

Clinical predictive modeling: Toward translational utility

The integrated four-gene molecular signature (apparent AUC = 0.867; nested cross-validated AUC = 0.68) provides a modest but hypothesis-generating model to classify stages based on molecular data. Consistent with this interpretation, a model combining the four-gene panel with the routine clinical variables (age, sex, smoking status, tumor location, differentiation grade) did not outperform the gene panel alone. This finding further suggests that the transcriptional panel, rather than these routine clinical variables, accounts for most of the stage-discriminative signal identified in this cohort. Identification of a stage-associated transcriptional threshold, including babA AUC 0.702, is consistent with - but does not establish - a broader concept of “microbial staging”. Under this concept, clinical assessment extends beyond detecting the pathogen to characterizing virulence transcript levels, which might eventually complement conventional staging[3]. A microbial “transcriptional score” could augment conventional TNM staging by providing a microbial “risk index”, enabling more rapid recognition of patients at high risk of disease progression even when morphological changes remain limited. We do not have evidence that this signature predicts rapid disease progression or clinical outcomes; the cross-sectional design captures association with stage at a single time point, not trajectory over time. However, this direction is aligned with current efforts to incorporate machine learning into personalized gastroenterology[27]. Before clinical interpretation is warranted, however, the transcriptional signature must be confirmed across diverse populations to determine whether it functions as a generalizable biological marker of GC.

Potential relevance to resource-limited settings

At the point of biopsy, the molecular classification framework could help identify patients with confirmed GC who are more likely to have advanced disease. This application should be tempered, however: Transcript-level qPCR requires RNA preservation, quality control, and molecular infrastructure that may not be readily available in the resource-limited settings where this tool could otherwise be most useful, including parts of the Pakistani healthcare system. Realizing this potential would require addressing that infrastructure gap alongside the validation steps noted above, before it could be considered practical for routine referral decisions.

Limitations and future directions

Interpretation of the findings requires consideration of several limitations: (1) Because the design is cross-sectional, causal relationships between bacterial transcriptional changes and malignant progression cannot be inferred; (2) Targeted qPCR was restricted to four virulence genes and therefore does not represent the full microbial or host transcriptomic landscape. Untargeted RNA sequencing or spatial transcriptomics in future studies may reveal additional regulatory complexity[28]; (3) Only 22 of the 200 cases overall were Stage III, including 18 of 159 in the H. pylori-positive sub-cohort, which is below the calculated power threshold of 44 participants per group (Supplementary material Section 3.3). Estimates for Stage III therefore have greater uncertainty than estimates for the other stages, and Stage IV included 37 of 159, both below the calculated power threshold of 44; (4) Specific canonical motifs, including the cagA EPIYA-C (Glu-Pro-Ile-Tyr-Ala) motif, have high background prevalence in this Pakistani cohort, consistent with our previous genomic work[8] - and may alter the predictive contribution of individual markers. Validation in geographically diverse populations is therefore necessary; (5) The analysis included only the 159 GC cases in which PCR amplified at least one H. pylori virulence gene and/or housekeeping gene. The 41 excluded cases differed significantly in histological differentiation (χ² = 8.77, P = 0.033), so the results should not be generalized to H. pylori-negative GC; (6) Host genetic and immune variables, including interleukin-1β polymorphisms and tumor necrosis factor-α levels, were not measured despite their established influence on H. pylori virulence gene expression and cancer risk. Unmeasured host factors may therefore confound part of the stage-associated signal; (7) The decision-tree classifier’s apparent AUC (0.867; sensitivity 74.5%, specificity 81.7%) came from fitting and evaluating on the same cohort; nested 5-fold cross-validation gave a lower, more defensible AUC of 0.68 (95%CI: 0.59-0.77); (8) The proportional-odds assumption held for cagA but not vacA, babA, or iceA (Brant test), so their odds ratios reflect an average effect across stage transitions rather than a uniform one; and (9) Of the 159 H. pylori-positive Tier two cases, 55 had non-adenocarcinoma gastric histologies, such as neuroendocrine tumors, lymphomas, and gastrointestinal stromal tumors. Because of this, TNM staging conventions may not apply equally to all cases in the pooled cohort. When the analysis was limited to the 104 cases with only adenocarcinoma, babA still showed a strong independent link to stage (OR = 0.583, 95%CI: 0.451-0.752, P < 0.001). In contrast, the association with iceA became weaker and was no longer statistically significant (OR = 0.837, 95%CI: 0.665-1.053, P = 0.129). This weaker association for iceA may be due to the smaller sample size and the resulting loss of statistical power in this restricted analysis. Nevertheless, when Stage III and Stage IV were combined in the sensitivity analysis, babA and iceA remained the same significant predictors with nearly unchanged coefficients, indicating that the stage-related signal discussed in limitation (3) is unlikely to arise solely from the small Stage III sample.

CONCLUSION

H. pylori virulence gene expression shows stage-associated variation in H. pylori-positive GC, with babA and iceA emerging as exploratory, stage-associated candidate biomarkers. The integrated four-gene signature outperformed any individual marker for distinguishing between stages within this cohort. Before this signature could be considered alongside histopathological staging, prospective validation in geographically diverse cohorts using this same normalization approach and direct comparison with established clinicopathologic staging models are needed.

ACKNOWLEDGEMENTS

We thank the patients and participating hospitals for contributing to this study. We also thank the laboratory and technical personnel at CUI for their support in sample processing and molecular assays.