An immune-related prognostic signature associated with immune landscape and therapeutic responses in gastric cancer

Data acquisition

The transcriptional expression data, somatic mutation patterns, and matching clinical information (including survival time, vital status, age, gender, tumor grade, and pathological stage) of GC patients were obtained from the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) database (Supplementary Table 1). An independent external GC cohort GSE84437 containing microarray data and corresponding clinical information was downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database as well (Supplementary Table 1). The gene expression profile and the detailed clinical annotations of an immunotherapy cohort (IMvigor210) were acquired online (http://research-pub.gene.com/IMvigor210CoreBiologies) and utilized to further validate the efficiency of the immune-related signature (Supplementary Table 2). The immune-related gene list was obtained from the ImmPort (https://www.immport.org) and InnateDB (https://www.innatedb.ca) databases (Supplementary Tables 3, 4). As all data were obtained from public databases, ethical approval was not necessary.

Data processing

For transcriptional expression data from the TCGA database, FPKM data were transformed into TPM (transcripts per kilobase of exon model per million mapped reads) values which is more similar to microarray data, to make them more comparable between different samples [10]. The mRNA Ensembl IDs in the TCGA dataset were switched to gene symbols according to GENCODE (https://www.gencodegenes.org/). For microarray data in GEO, the probe ID in the gene expression dataset was annotated to a gene symbol by platform files (GPL6947 Illumina HumanHT-12 V3.0 expression beadchip). Next, the above datasets were normalized by log2 transformation, and the batch effects were adjusted using the ComBat function in the R3.7.0 software “sva” package [11]. mRNAs with gene expression values of 0 in >90% samples both in TCGA and GEO datasets were excluded, as they were regarded as transcriptional noise. Averaging was performed for mRNAs with more than one-row expression values.

Differentially expressed gene analysis (DEGA), differentially expressed IRG (DE-IRG) screening and functional enrichment analysis

The R “limma” package was used for identifying differentially expressed genes (DEGs) between the adjacent normal and GC samples. Those with |logFC| > 1 and FDR < 0.05 were regarded as DEGs. Next, the DEGs and IRGs were intersected to obtain DE-IRGs, and a heatmap was used to visualize the gene expression profiles of DEGs and DE-IRGs respectively by R “pheatmap” package.

GO and KEGG enrichment analysis of DE-IRGs

Gene ontology (GO) and Kyoto Genome Encyclopedia (KEGG) enrichment analyses were conducted for the DE-IRGs via the R “ClusterProfiler” package to detect their underlying biological function. The GO terms and KEGG pathways with P < 0.05 were regarded as significant. The above results were visualized by R “ggplot2” and “enrichplot” packages eventually.

Weighted gene co-expression network analysis (WGCNA) of DE-IRGs

WGCNA is a method for exploring gene expression patterns of multiple samples [12]. Genes with similar expression patterns can be clustered, and the relationship between modules and specific traits or phenotypes can be analyzed. Firstly, Pearson’s correlation matrix was defined based on the interaction coefficients among genes. An adjacency matrix was defined with the threshold of Pearson coefficient exceeding 0.8, which was further used to construct a topological overlap matrix (TOM). Finally, the TOM matrix was applied for determining the co-expression gene modules, in this process modules with statistical significance (P < 0.05) were regarded as cancer-related modules, and genes involved in the modules were considered more important and used for subsequent analysis.

Identification of the immune-related signature (IRS)

In the discovery (TCGA) cohort, univariate Cox regression analysis was first conducted to detect the relationship between IRGs and GC prognosis. IRGs with P < 0.05 were regarded as having the potential to build the final prognostic signature. Next, to enhance the robustness of this final prognostic signature, the TCGA cohort was randomly divided into training and test subsets at a ratio of 5 to 5. In train subsets, IRGs with P < 0.05 in univariate Cox regression analysis were selected for inclusion in multivariate Cox regression analysis to establish the optimal prognostic signature. Particularly, the robustness of the prognostic signature was validated in the internal cohort (test cohort) and external cohort (validation cohort, GSE84437).

According to the prognostic signature, each patient was given a risk score by the following formula:

$$RiskScore={\displaystyle \sum _{k=1}^{n}Coe{f}_k\times B_k}$$

Where Coef_k is the coefficient and B_K is the normalized expression value of the IRGs included in the prognostic signature. GC patients were assigned to a low-risk group (LRG) or a high-risk group (HRG) according to the cohort-specific median risk score taken as the cut-off value. The performance of risk groups as determined by the risk score of the prognostic signature was evaluated by the Kaplan-Meier method with log-rank testing. Moreover, The ROC curves at 1, 3, and 5 years were plotted and AUC values were computed to estimate the predictive performance of the signature.

In addition, we divided patients into subgroups according to clinical characteristics such as age, pathological stage, etc. We then performed stratified survival analysis to study whether the IRS maintained predictive power in different cohorts.

Development and validation of a combined IRG-clinical nomogram

A nomogram is a statistical predictive tool that combines multiple prognostic factors to assess the survival probability for individual patients. To investigate whether the IRG prognostic signature possesses independent predictive capacity. First, risk scores together with clinical parameters including age, gender, tumor grade, and pathological stage were investigated by univariate Cox regression analysis to filter prognostic factors in the TCGA cohort. The above variables with P < 0.05 were then selected for inclusion in multivariate Cox regression analysis to determine independent prognostic ability. A combined prognostic model consisting of factors with P < 0.05 in multivariate Cox regression analysis was compared with the clinical model or IRS model respectively with regard to time-dependent AUC value and calibration curve. Finally, the model was visualized by nomogram and assessed by DCA to detect whether using this nomogram could yield closer associations with clinical net-benefit than other models.

Gene functions, tumor mutation burden, and somatic mutation profiling in the different risk subgroups

Gene set enrichment analysis (GSEA) is used to evaluate the distribution trend of genes in a predefined gene set in a gene table ranked by their relevance to phenotype, thereby judging their contribution to the phenotype [13]. The functional enrichment of genes in different risk subgroups was investigated in GSEA via the R “clusterProfiler” package. The somatic mutation characteristics of GC in the LRG and HRG were analyzed and visualized using the R “maftool” package. Moreover, the tumor mutation burden (TMB) for each patient was computed as mutations per million bases.

Correlations between different risk subgroups, TME and immune cell infiltration pattern

The ESTIMATE algorithm is a tool using expression data for the estimation of stromal cells and infiltrating immune cells in malignant tumors to predict tumor purity [14]. Therefore, utilizing ESTIMATE generates three scores: a stromal score (that captures the presence of stroma in tumor tissue), an immune score (that represents the infiltration of immune cells in tumor tissue), and the estimate score (that infers tumor purity) to evaluate the main cell types in the TME. Additionally, CIBERSORT, a deconvolution algorithm supported by R package with default parameters utilizing gene expression profiles to quantify immune infiltration, was used for evaluating the proportions of 22 immune cell types in each GC tumor tissue [15]. As CIBERSORT computes an empirical P-value for the deconvolution to denote the accuracy of results, we only retained those samples with a CIBERSORT P-value < 0.05 for subsequent analysis. Furthermore, to verify the results of CIBERSORT, other algorithms such as XCELL, TIMER, EPIC, and QUANTISEQ were used for analyzing the tumor infiltration immune cells as well [16–19].

Also, based on the single sample gene set enrichment analysis (ssGSEA) method supported by the R “GSVA” package, the 29 common immune-associated pathways for each sample were given an enrichment score to quantify this. Furthermore, immune cell abundance and their corresponding functional pathways were compared in different risk groups. Finally, survival analysis was performed for immune cells and immune functions respectively to comprehensively explore the relationship between immune infiltration and GC prognosis. The macrophages could play diverse roles in tumorigenesis and progression, such as M1 macrophages (prone to anti-tumor) and M2 macrophages (pro-tumor). Thus, the expression of markers for M0, M1, and M2 macrophages were analyzed in different risk groups. Besides, expression of immune checkpoint molecules PD-1, PD-L1, and CTLA4 in different risk groups were also analyzed. Moreover, due to the cytokine and chemokine are key factors for immune cell recruitment and functions, the expression of them in the risk model was analyzed as well.

The immunotherapy response and chemotherapeutic drug sensitivity in distinct risk subgroups

The tumor immune dysfunction and exclusion (TIDE) score calculated online (http://tide.dfci.harvard.edu/) was used to predict the immunotherapy response, the higher the TIDE score, the lower immunotherapy response [20]. Moreover, the IMvigor210 cohort was utilized to further verify the efficiency of the IRS in appraising immunotherapy responsiveness.

Based on the three public drug sensitivity databases, CGP (Cancer Genome Project), GDSC (Genomics of Drug Sensitivity in Cancer), and CTRP (Cancer Therapeutics Response Portal), the ridge regression models were built by the pRRophetic algorithm to predict chemotherapeutic responses. Gene expression profiling and risk grouping information were used in the model to estimate the half-maximal inhibitory concentrations (IC50) for each sample [21]. The prediction accuracy was evaluated by 10-fold cross-validation based on each training set. The smaller the IC50 value of the drug, the stronger its ability to inhibit cell growth and the better the effect of cancer treatment. Here, common GC therapeutic drugs, such as cisplatin, oxaliplatin, docetaxel, paclitaxel, 5-fluorouracil (5-Fu), capecitabine, and irinotecan were selected for analysis.

Immunohistochemical analysis

The protein expression data were acquired from the Human Protein Atlas (HPA) database, a largest and most comprehensive database for evaluating protein distribution in human tissues [22]. The protein expression of the IRGs signature in normal and GC tissues was determined using the immunohistochemical staining images, and then Image J was used to perform the quantitative analysis.

Validation of the immune-related signature by relative quantitative real-time PCR (qRT-PCR)

The expression levels of genes comprising the immune-relative signature were measured in a GC cell line (HGC-27, human gastric cancer cells) and a control cell line (GES-1, human gastric mucosal epithelial cells). All the cell lines were obtained from the National Infrastructure of Cell Line Resources (Beijing, China) and were in RPMI-1640 (FBS, Gibco, USA), 10% fetal bovine serum (FBS, Gibco, USA), and 1% penicillin/streptomycin (Gibco, Canada). All the cells were cultured at 37°C with 5% CO2. Total RNA was extracted from cells using the RNeasy Mini Kit (Qiagen, USA, Cat. 74104), and reverse transcription was subsequently performed using the 5∗All-in-one RT MasterMix (ABM, USA, Cat. No. G492). qRT-PCR was performed with a SYBR Green Real-time PCR Kit (Keygen Biotech, Nanjing, China, Cat. KGA1339-1) on a QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific, USA). All experiments were repeated at least three times. The RNA primer sequences are listed in Supplementary Table 5. Relative expression was calculated using the comparative threshold cycle (Ct) method. Simultaneously, the gene expression consisting of the signature was analyzed in the TCGA cohort.

Statistical analysis

All statistical analyses were conducted via R software (version 3.7.0). Continuous and categorical data were analyzed by Wilcoxon and Chi-square methods respectively. Survival was estimated by the Kaplan-Meier method and the statistical significance was determined by log-rank testing. Correlations between two continuous characteristics were analyzed by the Spearman method. Univariate and multivariate Cox regression analyses were performed using the R “survival” package. The time-dependent AUC value was calculated by the R “timeROC” package, and ROC curves were plotted by R “survivalROC” package. P < 0.05 was regarded as statistically significant.

Data availability statement

The data that support this study are openly available in online repositories and the raw data was uploaded (https://www.jianguoyun.com/c/sd/155e27d/6d60f576ace1d1cc).

Identification of DE-IRGs and functional enrichment analysis

With the cut-off set as logFC > 1, FDR < 0.05, there were 8833 DEGs identified totally, of which 1335 were down-regulated and 7498 were up-regulated (Supplementary Table 6). By searching ImmPort and InnateDB, 2483 and 1226 IRGs (Supplementary Tables 7, 8) were determined separately as well as 2660 IRGs were obtained by integrating IRGs from the two databases (Supplementary Table 9). Next, overlaps of DEGs and IRGs identified 493 DE-IRGs, comprising 184 that were down-regulated and 309 up-regulated (Supplementary Table 10). As shown in Figure 2A, 2B, there was a clear distinction of DEGs and DE-IRGs between normal and tumor samples.

Figure 2. Identification, functional enrichment, and WGCNA analysis of DE-IRGs. (A, B) Heatmap of DEGs and DE-IRGs, respectively. (C, D) GO enrichment and KEGG analysis for DE-IRGs. (E) The scale-free fit index for soft-thresholding powers. Left: the relationship between the soft-threshold and scale-free R2. Right: the relationship between the soft-threshold and mean connectivity. (F) Dendrogram and module colors for DE-IRGs. (G) The correlations heatmap between modules and sample types.

These 493 DE-IRGs were further analyzed via GO and KEGG to explore their function. In GO enrichment analysis, the biological process (BP) term showed that these DE-IRGs were enriched in leukocyte migration, regulation of the inflammatory response, cell chemotaxis, regulation of immune effector processes, and regulation of responses to external stimuli. The cellular component (CC) terms showed that these DE-IRGs were mainly enriched for the immunoglobulin complex, external side of plasma membrane, and collagen-containing extracellular matrix. Concerning molecular function (MF), these DE-IRGs were enriched in terms of receptor-ligand activity, cytokine activity, cytokine receptor binding, and growth factor activity (Figure 2C). Moreover, KEGG pathway analysis documented that these DE-IRGs were principally enriched in cytokine-cytokine receptor interaction, neuroactive ligand-receptor interaction, chemokine signaling pathway, viral protein interaction with cytokine and cytokine receptor, cAMP signaling pathway, and MAPK signaling pathways (Figure 2D).

Detecting hub DE-IRGs via WGCNA analysis

DE-IRGs were further selected for WGCNA analysis to distinguish different gene modules and then to identify the hub DE-IRGs. A soft threshold (power value) of 5 was selected to make sure that both the scale-free topology model fit index (R²) and mean connectivity were optimal (Figure 2E). A total of 5 modules was obtained by average linkage hierarchical clustering. Each module is shown in a different color in Figure 2F. Among these modules, MEturquoise possessed the highest correlation with GC traits (Figure 2G). There were 255 DE-IRGs included in the filtered MEturquoise module, which were chosen for further analysis (Supplementary Table 11).

Construction and verification of the immune-related prognostic signature

Utilizing univariate Cox regression analysis, a total of 36 genes was identified as being associated with GC prognosis (Figure 3A). We then constructed an 8-gene signature by multivariate Cox regression analysis in the training subset (Supplementary Table 12). The risk scores for each GC patients were calculated as follows:

Risk Score = (0.277) × RNASE2 + (0.255) × CGB5 + (0.501) × INGBE + (−0.348) × PTGER3 + (−0.341) × CTLA4 + (0.260) × DUSP1 + (0.071) × APOA1 + (0.234) × CD36.

Figure 3. Identification of the immune-related prognostic signature. (A) The forest plot of univariate Cox regression analysis for IRGs. (B–D) The Kaplan-Meier survival curves for high- and low-risk groups in train, test, and validation cohorts. (E–G) The ROC curve for the IRS in train, test, and validation cohorts.

Next, all GC patients were assigned to either a lower or higher risk group: LRG (n = 93) and HRG (n = 93). The results of Kaplan-Meier analysis and log-rank testing showed that patients in the LRG had a better overall survival rate (OS) than the HRG (Figure 3B). The AUC values of ROC curves for predicting survival outcomes at 1, 3, 5 years were 0.69, 0.86, and 0.87, respectively, demonstrating that the signature possesses good predictive capacity for GC prognosis (Figure 3C). In addition, the results of IRS validation showed that this signature maintain its prognosis predictive power in test subsets and validation cohorts as well (Figure 3D–3G).

Furthermore, we conducted stratified analysis for the IRS based on clinical characteristics such as age, gender, tumor grade, and pathological stage to explore the signature’s robustness. The results showed that patients in the HRG generally have a poor prognosis in different sub-cohorts compared with the LRG, which is consistent with the result from the integrated cohort. This confirms the robustness of this IRS in cohorts with different clinical features (Table 1).

Establishment and validation of an immune-clinical nomogram

The Cox regression analysis of the discovery cohort showed that age, pathological stage, and risk score were independent prognostic factors for GC patients (Table 2). Therefore, these factors were combined to construct a nomogram for predicting GC patients’ short- and long-term survival rates (Figure 4A). The time-dependent AUC value of the nomogram remained higher than for other models over time, suggesting that combining the IRS with age and pathological stage improved accuracy for predicting survival outcomes (Figure 4B). The calibration curve of the nomogram showed a good performance in consistency between prediction and actual observation, especially for 3-year OS. Also, the DCA curve suggested that the nomogram possessed better value for clinical applications than the other models (Figure 4C, 4D).

Figure 4. Construction of the immune-clinical nomogram. (A) The nomogram for predicting 1-year, 3-year, and 5-year OS for GC patients. (B, C) Time-dependent ROC curves for the nomogram, immune signature, age, and stage models at different time points in the TCGA and GEO datasets. (D, E) Calibration curves of observed and predicted probabilities for the nomogram in the TCGA and GEO datasets. (F, G) DCA curves for the nomogram in the TCGA and GEO datasets.

In the validation cohort, the time-independent AUC value and calibration curve of the nonogram model maintained its good performance for predicting patients’ OS as well (Figure 4E, 4F). DCA revealed that utilizing the nomogram could bring more clinical net benefit (Figure 4G). Thus, the nomogram comprising IRS and clinical characteristics (pathological stage and age) appeared highly accurate in predicting the short- and long-term OS of GC patients in both the discovery and validation cohorts. The above results demonstrated that using this immune-clinical nomogram for prognosis prediction assistance might result in significant clinical benefit.

GSEA analysis, TMB, and a gene mutation atlas in different risk subgroups

As shown in Figure 5A, 5B, GSEA analysis indicated that tumors in the HRG were enriched in complement and coagulation cascades, ECM receptor interaction, focal adhesion, neuroactive ligand-receptor interaction, and PPAR signaling pathway. In contrast, the LRG was enriched in aminoacyl tRNA biosynthesis, cell cycle, DNA replication, pyrimidine metabolism, and spliceosome terms.

Figure 5. GSEA, mutational landscape, and TMB in high- and low-risk groups stratified by the IRS. (A, B) The enriched pathways for different risk groups based on GSEA analysis (C, D) Waterfall plot of the top 20 mutant genes in the high- and low-risk groups. (E) Box plot for the TMB between the high- and low-risk groups. (F) Kaplan-Meier survival curve of high versus low TMB.

Summaries of gene mutation profiles for the different risk groups are shown in the Supplementary Figures 1 and 2, suggesting that the LRG has a higher overall mutation frequency than the HRG. Nonetheless, the top six genes with the highest mutation rates in the LRG were TTN (54%), TP53 (44%), MUC16 (35%), LRP1B (28%), ARID1A (28%), and SYNE1 (27%), which is similar to the HRG with TTN (40%), TP53 (39%), MUC16 (25%), LRP1B (20%), ARID1A (19%), and SYNE1 (16%). The most common mutation type was missense mutation in both LRG and HRG (Figure 5C, 5D).

TMB analysis in the HRG and LRG is shown in Figure 5E, 5F, suggesting that the latter has a higher TMB. In addition, survival analysis of patients stratified for TMB indicated that a high TMB resulted in a better prognosis than a low TMB.

Analysis of the tumor microenvironment, immune cell infiltration, and expression of markers for macrophages, immune checkpoint, cytokine, and chemokine in different risk groups

First, the difference between the TME and immune cell infiltration was analyzed for LRG and HRG. As shown in Figure 6A–6C, the stromal score and ESTIMATE score were both higher in HRG than LRG, whereas the immune score was not significantly different, indicating that the HRG had a higher proportion of stromal cells than LRG while tumor purity was lower.

Figure 6. TME and immune cell infiltration in different risk groups. (A–C) TME analysis based on ESTIMATE algorithm. From top to bottom: The stromal score, ESTIMATE score, and immune score. (D) Immune cell infiltration in high- and low-risk groups based on the CIBERSORT algorithm. (E) Immune-related pathways based on ssGSEA.

Second, by using the CIBERSORT, the abundance of immune cells in each GC sample was estimated and can be seen intuitively in the histogram (Supplementary Figure 3). The number of samples with significant levels of immune cell infiltration in the HRG was greater than in the LRG. Score difference analysis for immune cells’ abundance in different risk groups showed that compared with the HRG, the LRG had higher abundances of CD8+ T cells (P < 0.01^**), CD4+ activated memory T cells (P < 0.001^***), follicular T helper cells (P < 0.001^***), and M1 macrophages (P < 0.01^**). On the other hand, there was a lower abundance of monocytes (P < 0.001^***), M2 macrophages (P < 0.001^***), eosinophils (P < 0.05^*), and neutrophils (P < 0.001^***) (Figure 6D). Notably, the results of XCELL, TIMER, EPIC, and QUANTISEQ showed good consistency with the CIBERSORT (Supplementary Figure 4). Score differences for immune-related functional pathways suggested that the scores for cytolytic activity (P < 0.05^*), inflammation promoting (P < 0.01^**), MHC class I (P < 0.001^***), T cell co-inhibition (P < 0.05^*), Th1 cells (P < 0.05^*), and Th2 cells (P < 0.001^***) were higher in LRG than in HRG. In contrast, the scores for CCR (P < 0.01^**), DCs (P < 0.01^**), iDCs (P < 0.01^**), macrophages (P < 0.001^***), mast cells (P < 0.001^***), neutrophils (P < 0.001^***), T helper cells (P < 0.05^*), and type II IFN responses (P < 0.001^***) were lower in the LRG than in the HRG (Figure 6E). Moreover, survival analysis based on immune cells’ abundance in all samples implied that high infiltration of CD8+ T cells was connected with a favorable prognosis and that high abundance of M2 macrophage was related to an adverse prognosis (Figure 7A, 7B). This confirms that specific immune infiltration patterns have impact on patient prognosis. Survival analysis for immune-associated pathways suggested that high enrichment scores for cytolytic activity, inflammation promoting, T cell co-inhibition, and Th2 cell presence were associated with a favorable prognosis, whereas IDCs, mast cells, neutrophils, and type II IFN responses were associated with adverse prognosis (Figure 7C–7J).

Figure 7. Kaplan-Meier survival analysis for immune infiltration cell and immune-related pathway. (A, B) Kaplan-Meier survival curves for immune infiltrating CD8+ T cells and M2 macrophages. (C–J) Kaplan-Meier survival curves for immune-related pathway.

The expression analysis for markers of macrophages indicated that the marker gene of M0 macrophages like CD68 was highly expressed in the high-risk group (Supplementary Figure 5A), CD86 for M1 macrophages was lowly expressed in the high-risk group (Supplementary Figure 5B), similarly, NOS2 tended to be down-regulated in the high-risk group (Supplementary Figure 5C). Moreover, the markers for M2 macrophages like CD163 and CD206 were highly expressed in the high-risk group (Supplementary Figure 5D, 5E). The aforesaid results were consistent with the analysis of tumor infiltration immune cells for macrophages.

Analysis of immune checkpoints showed that PD-1, PD-L1, and CTLA4 were more highly expressed in LRG than in HRG, with decreasing expression as the risk score increased (Figure 8A–8F). The results of cytokine suggested that IFN-γ, IL-21, and OSM had higher expression in LRG than in HRG while IL-1α, IL-1β, IL-10, IL-24, TGF-β, EGF, and VEGF had lower expression in LRG than in HRG (Supplementary Figure 6A–6J). The analysis of chemokine showed that CCL2, CCL7, and CXCL8 were highly expressed in the HRG than in the LRG whereas the CXCL9, CXCL10, and CXCL11 were highly expressed in the LRG than in the HRG (Supplementary Figure 6K–6P).

Figure 8. The analysis of TIDE score and expression of immune checkpoints in high- and low-risk groups. (A–C) The expression of PD1, PD-L1, and CTLA4 in different risk groups. (D–F) The co-expression patterns between immune checkpoints and risk scores. (G–I) The scores of immune dysfunction, immune exclusion, and TIDE in different risk groups.

Immunotherapy response and chemotherapeutic drug sensitivity prediction based on the IRS to improve the GC patients’ survival

Compared with the HRG, the LRG presented with lower TIDE scores, indicating that the latter may have a greater response to ICI than the former (Figure 8G–8I). In the immunotherapy cohort, patients in the LRG had a markedly longer survival time (Figure 9A). Compared with HRG, the LRG had a better therapeutic advantage and immunotherapy response (Figure 9B). Moreover, the TMB was significantly elevated in LRG, which is closely linked to immunotherapeutic efficacy (Figure 9C). Also, the association between the IRS and survival on immunotherapy remained statistically significant after taking into account gender, smoking, ECOG score, immune phenotype, and TMB status (Figure 9D).

Figure 9. The IRS in the role of immune checkpoint blocker treatment. (A) Kaplan-Meier survival curve of the high- versus low-risk group in the immunotherapy cohort (IMvigor210 cohort). (B) The proportion of immune response to immunotherapy in high versus low-risk group. (C) The tumor mutation burden in the immunotherapy cohort was compared among distinct risk groups. (D) Multivariate Cox regression analysis of the IRS with features in the immunotherapy cohort.

Three drug response-related databases (CGP, GDSC, and CTRP) were utilized to investigate the association between the chemotherapeutic drug sensitivity and the IRS. Results suggested that patients in the LRG are generally more sensitive to chemotherapeutic drugs than those in the HRG (Figure 10).

Figure 10. Chemotherapy drugs’ sensitivity of the high- versus low-risk group. Differential analysis of IC50 for chemotherapy drugs in CGP (A–E), CTRP (F–I), and GDSC (J–P) databases.

Expression of the genes in the immune-related signature

The quantitative analysis of the immunohistochemical images determined that CD36, DUSP1, and PTGER3 showed lower protein expression levels in GC samples than in the adjacent normal tissues while CGB5 showed high levels in GC samples than in the normal tissues. The expression of APOA1, INHBE, and RNASE2 did not differ significantly in gastric cancer and normal samples, and the CTLA4 was not reported in the database (Supplementary Figures 7, 8).

Validation of the immune-related signature by qRT-PCR

The expression profiles of the eight genes comprising the prognostic signature were verified in GC and stomach cell lines by qRT-PCR. The result suggested that RNASE2, INHBE, CGB5, and CTLA4 were upregulated in GC cell lines, while PTGER2, DUSP1, CD36, and APOA1 were downregulated (Figure 11A), which showed good consistency with the expression analysis results in TCGA cohort (Figure 11B).

Figure 11. The expression of the immune-related signature in cell lines and TCGA cohort. (A) qRT-PCR results of the immune-related signature in GC cell lines (HGC-27) and control cell lines (GES-1). ^*P < 0.05, ^**P < 0.01, ^***P < 0.001. (B) The expression of the immune-related signature in TCGA cohort.