Characterization of Immune-Related Molecular Subtypes and a Prognostic Signature Correlating With the Response to Immunotherapy in Patients With Gastric Cancer

Abstract

Gastric cancer (GC) is a disease characterized by high molecular and phenotypic heterogeneity and represents a leading cause of cancer-related death worldwide. The tumor immune microenvironment (TIME) affects the response to immunotherapy and the prognosis of patients with GC. Explorations of the TIME in GC and characterization of molecular subtypes might enhance personalized treatment and facilitate clinical decision-making. In this study, two molecular subtypes were defined through unsupervised consensus clustering based on immune-related dysregulated genes. Then, patients with different molecular subtypes of GC were shown to have distinct differences in sensitivity to immune checkpoint blockers (ICBs). The immune-related prognostic signature was established utilizing least absolute shrinkage and selection operator (LASSO)-Cox regression analysis. Three independent external cohorts and the IMvigor210 cohort were introduced to validate the robustness of IPRS. scRNA-seq data of GC samples were used to decipher the underlying mechanisms of how IPRS contributes to the TIME. GC biospecimens were collected for RT-qPCR to further validate our findings. In summary, we characterized the abnormal TIME of GC and constructed a reliable immune-related prognostic signature correlating with the response to immunotherapy. This study may provide new strategies for developing individualized treatments for patients with GC.

Keywords: gastric cancer; gene signature; immunotherapy; molecular subtypes; prognosis.

Figure 1

The flow diagram of overall analysis.

Figure 2

Identification of IDGs and immune-related clusters. (A) Screening the IDGs by Venn diagram. (B, C) GO (B) and KEGG (C) enrichment analyses of IDGs. (D) The cumulative distribution function (CDF) curves for k = 2 to 6. (E) Three clusters were identified through consensus clustering. (F) PCA analysis revealed the dissimilarity among three clusters. BP, biological process; CC, cellular component; MF, molecular function.

Figure 3

Exploration of the TIME and characterization of immune-related molecular subtypes. (A, B) The fractions of tumor-infiltrating cells estimated by TIMER (A) and TIDE scores (B) among three clusters. (C, D) The fractions of tumor-infiltrating cells estimated by TIMER (C) and TIDE scores (D) between the two immune subtypes. (E) The levels of 28 immune cells of GC samples in the two immune subtypes were calculated through ssGSEA. (F) The relationships between immune subtypes and other previously reported subtypes in GC. NS. or ns, no statistical significance, *p < 0.05, **p < 0.01, ***p < 0.001.

Figure 4

Construction of an immune-related prognostic gene signature. (A) The changing trajectory of each variable in LASSO-Cox regression analysis. (B) Selection of the optimal lambda value with the minimum partial likelihood deviance. (C) Forest plot of multivariate Cox analysis of five selected genes based on the training set. (D) The coefficients of five selected genes. (E) Kaplan-Meier curves with log-rank p-values of the five selected genes. *p < 0.05, **p < 0.01.

Figure 5

Validation and assessment of the immune-related prognostic gene signature. (A–C) Kaplan-Meier curves of low and high IPRS groups in the training set (A), internal validation set (B), and the whole TCGA set (C). (D) ROC curves of the IPRS at 3 years and 5 years. (E) ROC curves of the IPRS and traditional clinical characteristics at 5-year-survival. (F) Time-dependent AUC curves of the IPRS and clinical factors.

Figure 6

The relationships between IPRS and clinicopathological features.

Figure 7

Investigation of the correlation between IPRS and tumor-infiltrating cells and the sensitivity to chemotherapies and ICBs. (A) IPRS was positively or negatively correlated with some tumor-infiltrating cell types. (B) The estimated IC50 for various chemotherapeutic drugs. (C) The distributions and proportions of potential responders and non-responders were estimated by TIDE scores between low and high IPRS groups. NS., no statistical significance, **p < 0.01, ***p < 0.001.

Figure 8

Exploration of the regulatory mechanisms resulting in the differences between low and high IPRS groups. (A, B) Hallmarks of cancers (A) and KEGG pathways (B) were investigated between low and high IPRS groups through GSEA. (C) IPRS was negatively correlated with the TMB. (D) IPRS showed significant discrepancies among different MSI statuses. MSI, Microsatellite Instable; MSS, Microsatellite Stable.

Figure 9

External validation of the immune-related prognostic gene signature. (A–C) Kaplan-Meier analyses of the low and high IPRS groups in three independent external cohorts on OS, including GSE14210 (A), GSE84437 (B), ACRG (C). (D–F) Kaplan-Meier analyses of the low and high IPRS groups in three cohorts on DFS, including TCGA-STAD (D), GSE14210 (E), ACRG (F). (G) Meta-analysis was conducted to evaluate the pooled HR of the immune-related prognostic signature. (H) IPRS in groups with a different anti-PD-L1 response status. (I) Kaplan-Meier analysis of patients with high and low IPRS in the IMvigor210 cohort. (J) The objective rate of clinical response to anti-PD-L1 immunotherapy in high and low IPRS groups in the IMvigor210 cohort. CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease. *p < 0.05.

Figure 10

Validation of the aberrant expression of the five selected IDGs in biospecimens. (A) The expression levels of five selected IDGs in adjacent tissues and gastric cancers obtained from TCGA. (B) RT-qPCR was performed using biospecimens to validate the expression pattern of five selected IDGs. *p < 0.05, **p < 0.01, ***p < 0.001.

Figure 11

Decipher the underlying mechanisms of how IPRS contributes to the TIME at the single-cell level. (A) UMAP plot showing sample origin after batch effect removal. (B) UMAP plot showing the expression levels of canonical marker genes for nine cell types. (C) Violin plot showing the expression distribution of canonical marker genes for nine cell types. (D) UMAP plot showing the distribution of nine cell types. (E) Bubble chart showing the expression patterns of genes involved in the prognostic model. (F) Bubble chart reflecting top ligand-receptor interactions of endothelial cells and endocrine cells communicating with other cell types. (G–I) Circle plots showing top ligand-receptor pairs in growth factor module (G), cytokine module (H), and immune checkpoint module (I) of cell-cell communication networks. Line width positively correlates with the expression levels of ligands; arrow width positively correlates with the expression levels of receptors.