Transcriptomic profiles and corresponding clinical information for STAD patients (n = 325) were downloaded from The Cancer Genome Atlas (TCGA) database. Additional validation cohorts were retrieved from the Gene Expression Omnibus (GEO), including GSE84437 (n = 423), GSE62254 (n = 297), GSE15459 (n = 175), and GSE26253 (n = 432). ComBat algorithm was used for batch effect correction using the “sva” package, followed by probe-to-gene mapping and z-score scaling to ensure comparability across datasets (Leek et al., 2012). In brief, we included histologically confirmed STAD primary tumors with available baseline transcriptomic data and overall survival (OS) information. We excluded samples lacking survival time or status, duplicate entries, and cases with OS follow-up <90 days (deaths and censorings, 8 cases excluded in the TCGA STAD cohort, 12 cases excluded in the GSE84437 cohort, 6 cases excluded in the GSE84437 cohort, 5 cases excluded in the GSE84437 cohort, and 8 cases excluded in the GSE84437 cohort) (Feng et al., 2023), to ensure a minimal, clinically meaningful follow-up window. Additionally, three immunotherapy datasets—GSE91061 (n = 98, melanoma), GSE78220 (n = 28, melanoma), and the IMvigor210 cohort (n = 298, urothelial cancer)—were utilized to evaluate the predictive value of the IES in determining the effectiveness of immunotherapy. Differentially expressed genes (DEGs) between STAD and normal tissues were identified using the “limma” package, with a threshold of |Log2FC| > 1.5 and p-value <0.05. A curated list of immune evasion-related genes (IRGs) was obtained from prior studies and public databases (Lawson et al., 2020; Wen et al., 2025), which are presented in Supplementary Table S1.

To screen for prognostic IRGs, univariate Cox proportional hazards regression was applied to the TCGA-STAD cohort. Genes significantly associated with overall survival (p < 0.05) were retained as candidate prognostic biomarkers for subsequent model construction. We then subjected the identified biomarkers to an integrative machine learning framework to develop a robust prognostic IES. This framework incorporated 10 machine learning techniques, including random survival forests, elastic net, Lasso, Ridge regression, stepwise Cox regression, CoxBoost, partial least squares regression for Cox models, supervised principal component analysis, generalized boosted regression modeling, and survival support vector machines. Following the framework outlined in previous studies (Liu et al., 2022; Li Z. et al., 2023), 101 algorithmic combinations were evaluated in the TCGA cohort via leave-one-out cross-validation. Each model was subsequently validated in GEO cohorts, and the Harrell’s concordance index (C-index) was calculated across all datasets. The model with the highest average C-index was selected as the optimal IES.

The optimal cut-off value for the IES score was determined using the “surv_cutpoint” function from the R package “survminer”. Patients were classified into high- and low-risk groups. Kaplan–Meier survival analysis and time-dependent ROC curves (via the “survivalROC” R package) were used to assess prognostic performance. We conducted both univariate and multivariate Cox regression analyses to identify prognostic factors in STAD patients. Based on these results, a prognostic nomogram was constructed using the “nomogramEx” R package, incorporating the IES-derived risk score and additional clinical variables. The predictive performance of the nomogram was assessed by comparing the predicted survival probabilities with the actual outcomes, which was illustrated using calibration plots.

To evaluate the immune microenvironment in STAD, we first applied the ESTIMATE algorithm to compute immune and stromal scores for each patient (Yoshihara et al., 2013). Additionally, seven other deconvolution tools—namely TIMER, xCell, MCP-counter, CIBERSORT, CIBERSORT-ABS, EPIC, and quanTIseq—were employed to estimate the composition of immune cells within tumor samples, selected for their complementary strengths (Li et al., 2020). The expression patterns of human leukocyte antigen (HLA)-associated genes and immune checkpoint regulators were illustrated using the “ggpubr” and “ggplot2” R packages. Furthermore, we conducted single-sample gene set enrichment analysis (ssGSEA) via the “GSVA” package to quantify immune-related functions and the infiltration levels of specific immune cell types.

To explore the predictive capability of the IES model in assessing immunotherapy response, we analyzed several key metrics: tumor immune dysfunction and exclusion (TIDE) (Fu et al., 2020), immunophenoscore (IPS) (Charoentong et al., 2017), and tumor mutation burden (TMB) (Samstein et al., 2019). TMB score of TCGA STAD patients were calculated with “maftools” package. Elevated TIDE and immune escape scores, alongside reduced IPS and TMB values, were considered indicative of stronger immune evasion and diminished likelihood of response to immune checkpoint blockade. Drug response for each STAD case was predicted based on half-maximal inhibitory concentration (IC50) values calculated using the “oncoPredict” R package (Maeser et al., 2021), referencing data from the Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013). Lower IC50 values represent increased sensitivity to therapeutic agents.

Protein expression levels of genes included in the IES were extracted from the Human Protein Atlas (https://www.proteinatlas.org/) (Colwill et al., 2011), which integrates proteomic and transcriptomic data to chart protein localization across human tissues and tumor types. Expression profiles were obtained from both normal tissue and cancer-specific datasets.

Normal gastric epithelial (RGM-1) and STAD cell lines (NCI-N87, CRL-5822, RPMI1640, BGC-823, and HGC-27) were acquired from the Shanghai Institute of Biochemistry and Cell Biology. All cells were maintained in ATCC-recommended media supplemented with fetal bovine serum (FBS; Gibco) and 1% penicillin-streptomycin (Sigma-Aldrich), and cultured at 37 °C in a humidified atmosphere containing 5% CO₂. For gene silencing experiments, NCI-N87 and HGC-27 cells were transfected with KLF16-targeted siRNA or negative control siRNA using Lipofectamine 3000 (Invitrogen) according to the manufacturer’s instructions.

Total RNA was isolated using TRIzol reagent (Takara Bio), followed by cDNA synthesis using oligo(dT) primers. Quantitative real-time PCR (RT-qPCR) was performed on an ABI 7900HT system (Thermo Fisher Scientific) with SYBR Premix Ex Taq (Takara Bio), and gene expression levels were normalized to GAPDH.

For proliferation assessment, gastric cancer cells were seeded into 96-well plates at a density of 5,000 cells per well. At indicated time points, Cell Counting Kit-8 (CCK-8; Beyotime) reagent was added, and optical density (OD) values were measured. The proliferation index was calculated as the OD ratio at each time point relative to baseline.

A total of 2962 DEGs between STAD and normal tissues were identified (Supplementary Figure S1A). By overlapping these DEGs with previously defined IRGs, we obtained a subset of 40 candidate genes (Supplementary Figure S1B). Univariate Cox regression analysis further revealed 8 IRGs (KLF16, ANXA5, TAP1, DOT1L, PRKCSH, NRP1, MARCKS, and AKR1B1) significantly associated with overall survival in STAD, which were considered as potential prognostic biomarkers (Supplementary Figure S1C).

To construct a robust prognostic model, we applied ten different machine learning algorithms and developed 101 predictive models across TCGA and GEO cohorts. The model generated by the LASSO algorithm achieved the highest average C-index (C-index = 0.82) and was selected as the optimal IES (Figure 1A). In LASSO algorithm, four genes (KLF16, ANXA5, TAP1, and DOT1L) were selected from above 8 potential prognostic biomarkers for the construction of optimal IES. The IES score (risk score) was calculated with the following formula: IES score = (0.282) × KLF16^exp + (−0.098) × ANXA5^exp + (−0.125) × TAP1^exp + (−0.035) × DOT1L^exp. Based on the optimal cut-off value (0.035), STAD patients were divided into high- and low-risk groups. Kaplan–Meier survival analysis demonstrated significantly poorer overall survival in the high-risk group across all validation cohorts (TCGA, GSE84437, GSE62254, GSE15459, and GSE26253) (Figures 1B–F). The corresponding ROC curves indicated that the IES exhibited excellent predictive accuracy for 1-, 3-, and 5-year survival (Figures 1G–K).

We compared the prognostic performance of the IES with traditional clinical factors such as age, gender, and stage. The IES consistently demonstrated a higher C-index across all datasets (Figure 2A). Both univariate and multivariate Cox regression analyses confirmed the IES as an independent risk factor for STAD prognosis (Figures 2B,C). A prognostic nomogram incorporating the IES and clinical parameters was constructed (Figure 2D), and calibration plots revealed excellent concordance between predicted and observed survival outcomes (Figure 2E).

To investigate the immunological implications of IES, we assessed the immune cell infiltration profiles using seven deconvolution algorithms. The IES-based risk score was negatively correlated with immune infiltration, including CD8⁺ T cells and dendritic cells, while positively correlated with immunosuppressive M2 macrophages (Figures 3A–D). Patients in the low-risk group exhibited higher enrichment scores for immune cells and immune-related functions, such as iDCs, mast cells, NK cells, TILs, cytolytic activity and T cell co-stimulation (Figures 3E,F). Moreover, they had significantly higher immune, stromal, and ESTIMATE scores, indicating a more active immune microenvironment (Figure 3G).

We further explored the role of IES in predicting immunotherapy benefits. Patients in the low IES group showed increased expression of immune checkpoint molecules and HLA genes (Figures 4A,B), elevated PD1 & CTLA4 immunophenoscore (Figure 4C), and higher TMB scores (Figure 4D). Elevated expression of HLA-related genes and TMB score predicted a higher chance of immunotherapy benefits (Lin and Yan, 2021; Hodi et al., 2021). A stronger response to immunotherapy was indicated by low TIDE and immune escape scores (Fu et al., 2020; Lin et al., 2020). Conversely, these patients had significantly lower TIDE scores, immune escape scores, and intra-tumor heterogeneity (ITH) scores (Figures 4E–G), suggesting a reduced potential for immune evasion. In three independent immunotherapy cohorts (GSE91061, GSE78220, IMvigor210), lower IES scores were consistently associated with better response rates and improved survival outcomes (Figures 4H–J).

We evaluated the association between IES and drug sensitivity. Patients with low IES-based risk score exhibited significantly lower IC50 values for multiple chemotherapeutic agents (e.g., Cisplatin, Oxaliplatin, Docetaxel, Gemcitabine) and targeted therapies (e.g., Lapatinib, Erlotinib, Dasatinib, Afatinib), suggesting greater sensitivity to conventional treatments (Figures 5A,B).

To gain mechanistic insights, we performed functional enrichment analyses. High-risk patients showed enhanced activation of hallmark cancer pathways, including mTORC1 signaling, hypoxia, NOTCH signaling, glycolysis, angiogenesis, IL2-STAT5 signaling, G2M_checkpoint, and hedgehog signaling, indicating a more aggressive tumor phenotype (Figure 6).

Among the key genes comprising the IES, KLF16 was selected for further validation. Immunohistochemistry confirmed overexpression of KLF16 in STAD tissues (Figure 7A), and RT-qPCR showed elevated mRNA levels of KLF16 in most of STAD cell lines compared to normal gastric epithelial cells (Figure 7B). Knockdown of KLF16 significantly inhibited the proliferation of STAD cells, supporting its oncogenic role (Figure 7C).