Identification of immunotherapy-related subtypes, characterization of tumor microenvironment infiltration, and development of a prognostic signature in gastric carcinoma

Data acquisition

The immunotherapy cohort of gastric carcinoma (PRJEB25780) was obtained from the Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) database to identify immunotherapy-related genes [16, 17]. The RNA-sequencing profiles and corresponding clinicopathological information of GC were obtained from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) [18]. Detailed information of somatic mutation and copy number variation (CNV) data files were retrieved for further analysis. The immunophenoscore (IPS) of GC samples was obtained from the Stomach Adenocarcinoma (STAD) project of The Cancer Immunome Atlas (TCIA, https://tcia.at/) [19]. Four cohorts (GSE15459 [20], GSE84437 [21], GSE62254 [22], and GSE26253 [23]) were collected from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) to validate the predictive efficiency of the risk score signature in prognosis prediction [24]. IMvigor210 cohort includes RNA-sequencing profiles and detailed clinical parameters of urothelial carcinoma patients who received anti-PD-L1 therapy [25]. GSE78220 is an immunotherapy cohort where patients with melanoma were treated with anti-PD-1 agents [26]. We collected IMvigor210 and GSE78220 cohorts for evaluating the predictive efficiency of the risk score signature in immunotherapy response. The raw count data obtained from GEO database were normalized by applying the “limma” package in R [27].

Identification of IRGs in GC

The IRGs were defined as genes that were differentially expressed between immunotherapeutic responder and non-responder groups. In our study, patients with GC in PRJEB25780 cohort were divided into responder and non-responder groups according to their responsiveness to anti-PD-1 immunotherapy. Then, the “limma” package in R was applied to screen out IRGs, the cutoff conditions were set to adjusted P-value < 0.05 and absolute value of log2 fold change (log2 FC) ≥ 0.585. Besides, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis to further clarify the potential functional annotation of the IRGs based on the Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/) database [28]. Somatic mutations and CNV are important events leading to genetic alterations. Thus, we first analyzed the overall genetic alteration atlas of the IRGs based on cBioPortal (https://www.cbioportal.org) platform [29]. Then, the detailed mutation and CNV information of the IRGs were analyzed.

Consensus clustering analysis of IRGs

The “ConsensusClusterPlus” package with 1000 repetitions was used for performing unsupervised clustering analysis to identify distinct molecular subtypes in the TCGA cohort based on IRGs expression [30]. Then, differences in prognosis and clinicopathological parameters between subgroups were analyzed.

Correlations between the subtypes and TME infiltration

The “CIBERSORT” algorithm was used to quantify the infiltrated abundances of 22 immunocytes for each GC patient with TCGA expression data [31]. The “ESTIMATE” algorithm was used to evaluate the purity of tumors from different clusters [32]. We compared the differences in immunocyte and immune score between different clusters to verify the immune characteristics of the IRG clusters. In addition, we applied a series of indexes which was usually used for predicting the response to immunotherapy, including tumor mutation burden (TMB), immune checkpoint biomarkers (ICBs), TIDE score, microsatellite instability and IPS, to assess the correlation between the subtypes and the effect of immunotherapy.

Creating and confirming the predictive risk score signature

Differentially expressed genes (DEGs) between the immunotherapy-related subgroups were identified using the “limma” package in R, the threshold was set as adjusted P-value < 0.05 and absolute value of log2 FC ≥ 0.585. Similarly, gene ontology (GO) enrichment analysis and KEGG pathway enrichment analysis were conducted to clarify the pathways that were considerably enriched. Subsequently, univariate Cox and LASSO regression analyses were carried out to screen out the optimal prognostic biomarkers among these DEGs and included them in the immunotherapy-related risk score (IRRS) signature. The following formula can calculate the risk score of all cases: $RiskScore={\displaystyle \sum _1^n\text{Exp}(i)*Coef(i)}$, where Exp is the expression value of the gene, and Coef is the LASSO regression analysis coefficient of each gene in the signature. Kaplan-Meier survival analysis was utilized to explore the predictive ability of the IRRS in GC prognosis. Area under the receiver operating characteristic (ROC) curve was utilized to evaluate the diagnostic efficacies. Univariate and multivariate Cox proportional hazards regression analyses were used to assess whether the signature could be served as an independent prognostic factor. Besides, the correlations between the IRRS and the clinicopathological parameters, including age, gender, grade, stage, and microsatellite instability, were determined using chi-square tests.

Importantly, the risk score was also calculated in four external cohorts (GSE15459, GSE84437, GSE62254, and GSE26253) to validate the predictive ability of the IRRS signature in distinct populations.

Establishment of a nomogram

To enhance the clinical utility of the risk signature, a nomogram containing the IRRS and other independent prognostic factors was constructed using the “rms” package in R, and then was applied to predict the 1-, 3-, and 5-year survival rate of GC patients. ROC and calibration curves were used to evaluate the accuracy of the nomogram.

The role of IRRS in the prediction of immunotherapeutic benefits

We compared the differences in several immunotherapeutic response indexes (ICBs expression, TMB score, TIDE score, microsatellite instability and IPS) between the risk groups to assess the predictive ability of the risk signature in the prediction of immunotherapeutic benefits. Patients in IMvigor210 and GSE78220 cohorts were classified into responder (including partial response (PR) and complete response (CR)) and non-responder (including progressive disease (PD) and stable disease (SD)) groups according to the patients’ response to immunotherapy. Then, we calculated the risk score of each patient in IMvigor210 and GSE78220 cohorts based on the formula generated in TCGA cohort and analyzed its impact on the prognosis and the efficacy of immunotherapy.

Statistical analysis

All visualization and statistical analyses were performed by using R software (version 4.2.1, https://www.r-project.org/) and the corresponding feature packages. Two-sided P-value < 0.05 was considered as significant thresholds for all statistical tests.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/), and Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) projects.

Identification of IRGs in GC

The RNA-sequencing profiles and corresponding clinical parameters of GC patients in PRJEB25780 dataset were obtained from the TIDE database. Then, a total of 63 genes that differentially expressed between the responder and non-responder groups were identified as IRGs (Supplementary Table 1). The volcano plot of the IRGs was displayed in Figure 1A, 1B presented the expression heatmap of these IRGs. KEGG pathway annotation analysis revealed that the IRGs were primarily enriched in Pathways in cancer, Cytokine-cytokine receptor interaction, Rap1 signaling pathway, and PD-L1 expression and PD-1 checkpoint pathway in cancer (Figure 1C).

Figure 1. Identification of IRGs. (A) Volcano plot of IRGs. The red dots represent the upregulated genes, the green dots represent the downregulated genes, and the black dots represent genes with no significant difference in expression. (B) Expression heatmap of the IRGs. Red represents upregulated genes, and blue represents downregulated genes. (C) KEGG enrichment analysis of the IRGs.

Genetic variation landscapes of IRGs in GC

The somatic mutation and CNV frequencies of the 63 IRGs in GC patients were analyzed based on TCGA cohort. As presented in Figure 2A, 2B, genetic alterations of the IRGs occurred in 302 of 434 (69.59%) GC patients. Among them, FAT4 (19%), DCHS1 (7%), NID1 (5%), and JAKMIP1 (5%) possess the highest mutation frequency, and their main mutation types are missense mutation, frame shift del, and nonsense mutation. The changes in IRGs with CNV features on chromosomes was showed in Figure 2C, 57 out of the 63 IRGs had frequent copy number alterations, and most of the IRGs were accumulated on copy number loss rather than copy number gain (Figure 2D).

Figure 2. The alterations of IRGs in GC. (A) Gene alteration frequency of IRGs in GC based on TCGA project. (B) Landscape of genomic aberrations of the top 20 IRGs in GC. (C) The sites of CNV variation in IRGs on the 23 chromosomes. (D) Frequencies of CNV gain, loss, and non-CNV among IRGs.

Identification of IRG subtypes in GC

An unsupervised clustering strategy was used to classify 371 GC samples into cluster 1 (n=184) and cluster 2 (n=187) subgroups based on the expression patterns of 63 IRGs (Figure 3A). The principal component analysis (PCA) revealed that the discrimination between cluster 1 and cluster 2 is great (Figure 3B). Kaplan-Meier survival analysis revealed that patients in cluster 1 possess better OS and progression free survival (PFS) rate than those in cluster 2 (Figure 3C, 3D). The correlations between the clinical features and expression patterns of the IRGs was presented in the heatmap (Figure 3E). Patients of cluster 1 had a higher proportion of microsatellite instability-high (MSI-H) and alive outcome than cluster 2.

Figure 3. Identification of IRG subtypes in GC. (A) Consensus matrix heatmap defining two clusters (k = 2) and their correlation area. (B) A considerable transcriptome divergence between the two subtypes is seen by PCA analysis. (C) The Kaplan–Meier curve of OS analysis between different cluster groups. (D) The Kaplan–Meier curve of PFS analysis between different cluster groups. (E) Differences in clinicopathologic features and expression levels of IRGs between the two subtypes.

Analysis of TME infiltration in distinct subtypes

We carried out the “CIBERSORT” and “ESTIMATE” methods to determine the difference in TME infiltration, including 22 subtypes of immune cells and TME scores, between the two clusters in order to learn more about how IRGs work in the TME. As presented in Figure 4A, the difference in TME infiltration between these two clusters is great. Samples in cluster 1 seemed to exhibit remarkably lower stromal scores, lower immune scores, and higher tumor purity scores, compared with those of cluster 2 (Figure 4B). Kaplan-Meier survival analysis revealed that high stromal scores and low tumor purity scores were remarkably associated with poor prognosis in GC patients (Figure 4C). With regard to immune cell infiltration, the cluster 1 subgroup was characterized by the high infiltration of follicular helper T cells, activated memory CD4+ T cells, CD8+ T cells and M1 macrophages, whereas the cluster 2 was characterized by the high infiltration of resting memory CD4+ T cells, naive B cells, resting Dendritic cells, Monocytes, resting Mast cells, and Eosinophils (Figure 4D, 4E). Among them, elevated infiltration abundances of CD8+ T cells, follicular helper T cells and activated memory CD4+ T cells were significantly correlated with favorable clinical outcome of GC, while elevated infiltration abundance of naive B cells was unfavorable (Figure 4F).

Figure 4. Immune landscape of IRG subtypes in the TCGA cohort. (A) Differences in TME scores and infiltration levels of immune cells between the two distinct subgroups. (B) The comparisons of stromal score, immune score, and tumor purity between different subgroups. (C) The Kaplan–Meier survival regarding the stromal score, immune score, and tumor purity in GC patients. (D) The boxplot of the differences of the immune cells infiltration between two distinct subgroups. (E) Relative proportion of immune infiltration in two distinct subgroups. (F) The Kaplan–Meier survival analysis of the correlation between the infiltration levels of the immune cells and OS in GC. *P < 0.05, **P < 0.01, ***P < 0.001.

Subsequently, we determined whether the IRG subtypes had a significant correlation with the immunotherapy effect. Patients in cluster 1 possess higher TMB score than those in cluster 2, whereas the cluster 2 possess higher TIDE score (Figure 5A, 5B). Besides, the cluster 1 was characterized by the high proportion of MSI-H and low proportion of microsatellite stability (MSS), compared with cluster 2 (Figure 5C). In terms of ICBs, CD274 (also known as PD-L1) and LGALS9 expression in cluster 1 were significantly elevated, while the expression of CD276, HAVCR2, TIGIT, and PDCD1LG2 in cluster1 were downregulated (Figure 5D). Moreover, we calculated the IPS (CTLA4-/PD-1-, CTLA4+/PD-1-, CTLA4-/PD-1+ and CTLA4+/PD-1+) to reveal the immunogenicity of patients in each subgroup. As a result, the IPS were remarkably elevated in the cluster 1 subgroup, which appeared to have stronger immunogenicity. (Figure 5E). In sum, these findings implied that effective immune response was more activated in the cluster 1 group, which mean that the cluster 1 group exhibited a better response to immunotherapy.

Figure 5. The estimation of two IRGs subtypes in immunotherapy response. (A) The difference in TMB between two distinct subtypes. (B) The difference in TIDE score between two distinct subtypes. (C) Relative proportion of microsatellite instability in two distinct subtypes. (D) The expression of ICBs in two distinct subtypes. (E) The difference of IPS in two distinct subtypes.

Construction of the IRRS signature

As listed in Supplementary Table 2, 3065 DEGs (|log FC| ≥ 0.585, adjusted P < 0.05) between cluster 1 and cluster 2 subgroups were identified, and the expression heatmap of these DEGs was presented in Figure 6A. GO annotation analysis revealed that the DEGs are significantly annotated in immune-associated crosstalk, such as immune response, adaptive immune response, B cell receptor signaling pathway, and phagocytosis (Figure 6B and Supplementary Table 3). KEGG pathway annotation analysis found that the DEGs are primarily enriched in several areas of carcinogenesis-associated pathways, including ECM-receptor interaction, Focal adhesion, PI3K-Akt signaling pathway, and Pathways in cancer (Figure 6B and Supplementary Table 4). Subsequently, 625 genes presented significant correlations with the OS of GC patients were filtered out for further research using univariate Cox analysis (Supplementary Table 5), and three of which were eventually screened out for constructing the IRRS signature via utilizing LASSO regression analysis with minimized lambda (Figure 6C, 6D). The forest plot illustrated the correlations between the expression levels of the three genes in the IRRS signature and the prognosis of GC patients (Figure 6E). As presented in Figure 6F, we applied the following equation to calculate the risk score of each patient: IRRS = (Exp_NRP1 × 0.0354) + (Exp_SERPINE1 × 0.0478) + (Exp_CPNE8 × 0.0749). According to the median value of the risk score, 371 samples in the TCGA GC cohort were divided into the low- (n = 186) and high-risk (n = 185) groups, and the association of IRRS, IRG subtypes, and vital status of these patients was showed in the Sankey plot (Figure 6G).

Figure 6. Construction of the IRRS model in the TCGA cohort. (A) Expression heatmap of the DEGs between the two distinct subtypes. (B) GO and KEGG enrichment studies of DEGs between different cluster groups. (C) LASSO coefficient profiles of 3 selected genes in the 10-fold cross-validation. (D) Partial likelihood deviance was revealed by the LASSO regression model in the 10-fold cross-validation. (E) Forest plot of hazard ratios for three selected prognostic variables. (F) The coefficients of the three prognostic variables in the LASSO regression model. (G) Sankey plot displayed the correlations among the cluster subgroups, IRRS and vital status.

IRRS for the prognostic prediction of GC

We determined the clinical significance of the IRRS in GC based on TCGA cohort. The Kaplan–Meier survival analysis revealed that patients in the high-risk subgroup possessed a significant poorer OS than patients in the low-risk subgroup (Figure 7A). The ROC curve for 5-year OS of the IRRS signature exhibited an area under ROC (AUC) value of 0.721, implying that the signature had moderate sensitivity and specificity (Figure 7B). Patients in the high-risk subgroup possessed a higher proportion of histologic G3, MSS, and deceased outcome than those in the low-risk subgroup, and the results of the chi-square test revealed that NRP1, CPNE8 and SERPINE1 were all overexpressed in the high-risk subgroup (Figure 7C). Besides, we found that patients’ risk of death increased with increasing risk score (Figure 7D). Furthermore, we conducted a subgroup analysis on the basis of age, gender, histologic grade, and TNM stage. As expected, the OS of GC patients was better in the low-risk subgroup than in the high-risk subgroup in each subgroup (Figure 7E).

Figure 7. Clinical significance of the IRRS model in TCGA cohort. (A) The Kaplan–Meier survival analysis of the signature for predicting the OS of patients in TCGA cohort. (B) Time-dependent ROC analysis of the signature for predicting the OS of patients in TCGA cohort. (C) Differences in clinicopathologic features and expression levels of prognostic variables between the low- and high-risk groups. (D) The distribution of the IRRS and the vital status of patients in the TCGA cohort. (E) The subgroup survival analysis according to the age, gender, histologic grade, and tumor stage.

Validation of the IRRS in external cohorts

We performed validation analysis for evaluating whether the IRRS has clinical application value in distinct populations. First, the risk score for individual patients in GSE15459, GSE84437, and GSE62254 cohorts was calculated according to the calculation formula derived from TCGA cohort. Similarly, GC patients from these validation cohorts were divided into low- and high-risk subgroups according to the median value of the risk score (Figure 8A–8C). As shown in Figure 8D–8F, obvious differences in the OS probability were observed between the low- and high-risk subgroups in each validation cohort, GC patients with higher risk score seemed to exhibit poorer prognosis. In addition, we further determined whether the IRRS has predictive ability in relapse-free survival (RFS) of GC patients from GSE26253 cohort. Consistent with the results in OS, patients with higher risk score possessed higher rates of tumor recurrence (Figure 9). To sum up, these findings indicated the reliable capability of the IRRS in distinguishing differences in clinical outcomes of GC.

Figure 8. External validation of the IRRS model in predicting OS of GC patients based on three independent cohorts. (A) The distribution of the risk score, the vital status of patients, and the expression heatmap of three prognostic variables in the GSE15459 cohort. (B) The distribution of the risk score, the vital status of patients, and the expression heatmap of three prognostic variables in the GSE84437 cohort. (C) The distribution of the risk score, the vital status of patients, and the expression heatmap of three prognostic variables in the GSE62254 cohort. (D) The Kaplan–Meier survival analysis of the IRRS signature for predicting the OS of patients in the GSE15459 cohort. (E) The Kaplan–Meier survival analysis of the IRRS signature for predicting the OS of patients in the GSE84437 cohort. (F) The Kaplan–Meier survival analysis of the IRRS signature for predicting the OS of patients in the GSE62254 cohort.

Figure 9. External validation of the IRRS model in predicting RFS of GC based on the GSE26253 cohort. (A) The Kaplan–Meier survival analysis of the IRRS signature for predicting the RFS of GC patients based on the GSE26253 cohort. (B) The distribution of the risk score of GC patients in the GSE26253 cohort. (C) The expression heatmap of three prognostic variables in the GSE26253 cohort. (D) The distribution of the recurrent status of GC patients in the GSE26253 cohort.

The IRRS can be served as an independent risk predictor in GC

To investigate whether the IRRS can be served as a clinically independent risk predictor for GC patients, the IRRS, age, gender, histologic grade, and TNM stage were enrolled as covariates to perform the univariate and multivariate Cox analyses. The results indicated that age, TNM stage and IRRS are independent prognostic predictor that could be used for predicting the prognosis of GC patients (Figure 10A, 10B). Then, a nomogram was built to improve the clinical power of the IRRS by combining the prognostic factors above (Figure 10C). Multi-variables ROC analysis indicated that the nomogram had optimum predictive performance compared with other single factors (Figure 10D). Time-dependent ROC of 1-, 3-, and 5-year OS also revealed that the nomogram had high predictive power (Figure 10E). Besides, calibration charts revealed that the predictive capability of the nomogram in 1-, 3-, and 5-year periods was highly accurate, confirming its utility in forecasting the clinical outcome of GC patients (Figure 10F).

Figure 10. Nomogram developed for predicting the probability of 1-, 3- and 5-year OS in TCGA cohort. (A) Univariate Cox analysis containing IRRS and clinicopathological parameters. (B) Multivariate Cox analysis containing IRRS and clinicopathological parameters. (C) The comprehensive nomogram for predicting probabilities of GC patients with 1-, 3- and 5-year OS in TCGA dataset. (D) A comparison of ROC curve showed the superiority of the nomogram. (E) ROC curves chart of the comprehensive nomogram predicting the 1-, 3- and 5-year survival rate. (F) Calibration plot of nomogram for predicting probabilities of 1-, 3-, and 5-year survival of GC patients. Nomogram-predicted probability of survival is plotted on the x-axis; actual survival is plotted on the y-axis.

Risk score-based treatment strategy for GC

The expression of four ICBs (PD-L1, HAVCR2, TIGIT, and CTLA4) were compared between the low- and high-risk subgroups. As shown in Figure 11A, the results found that the presented immunomodulators were significantly upregulated in the high-risk subgroup. Subsequently, we determined the association between the IRRS and IPS (CTLA4-/PD-1-, CTLA4+/PD-1-, CTLA4-/PD-1+ and CTLA4+/PD-1+), which are recognized indexes to predict patients’ responses to ICIs by evaluating the immunogenicity. The results demonstrated that the IPS were remarkably elevated in the low-risk subgroup, implying that GC patients in the low-risk subgroup might be more sensitive to anti-cancer immunotherapy (Figure 11B). We also compared the distribution of TMB and TIDE scores between the two IRRS subgroups and revealed that the TMB score was lower in the high-risk subgroup, and the TIDE score was contrary, indicating that patients in the high-risk subgroup might respond worse to anti-cancer immunotherapy (Figure 11C, 11D). Besides, we explored the correlation between IRRS and microsatellite instability, and revealed that patients in the low-risk subgroup possess a higher percentage of MSI-H, and MSS more happened in patients in the high-risk subgroup (Figure 11E). In sum, all these results implied that GC patients in low-risk subgroup exhibited a better response to anti-cancer immunotherapy, IRRS might be a potential biomarker for immunotherapy of GC.

Figure 11. The immunotherapeutic benefit of the IRRS model. (A) The expression of four ICBs (PD-L1, HAVCR2, TIGIT, and CTLA4) in low- and high-risk groups. (B) The difference in IPS scores between low- and high-risk groups. (C) Box plots and scatter diagram of the TMB score and IRRS. (D) Box plots and scatter diagram of the TIDE score and IRRS. (E) Boxplot and Bar diagram of the microsatellite instability and IRRS. (F) Kaplan–Meier curve of OS for patients with high and low IRRS subtypes in IMvigor210 cohort. (G) Boxplot and Bar diagram displayed the response to immunotherapy in low- and high-risk groups in IMvigor210 cohort. (H) Boxplot and Bar diagram displayed the response to immunotherapy in low- and high-risk groups in GSE78220 cohort.

Subsequently, IMvigor210 and GSE78220 cohorts were applied as external validation cohort for verifying the predictive ability of the IRRS in anti-tumor immunotherapy. In the IMvigor210 cohort, patients in the low-risk subgroup had remarkably better prognosis than those in the high-risk subgroup (Figure 11F).

Besides, we analyzed the differences in immunotherapy response subtypes (CR, PR, SD, and PD) between different IRRS subgroups and found that patients with lower risk score possess a higher percentage of CR/PR, and SD/PD more happened in patients with higher risk score (Figure 11G, 11H). These findings further confirmed that the IRRS can be served as a great biomarker to forecast the sensitive to immunotherapy.