The landscape and prognostic value of tumor-infiltrating immune cells in gastric cancer

Data profile

The gene expression profiles were obtained from TCGA and ICGC datasets (downloaded in April 2019), and then subjected to background correction and normalization with Perl 5.0 (http://www.perl.org/). Meanwhile, relevant clinical characteristics of cancer cases were also collected. Patients with a follow-up time <30 days or a lack of pathological diagnosis were excluded from the study. The differences of demographic and baseline characteristics between two sets were compared by 2 tests, using SPSS 20.0 (IBM Corp.). P < 0.05 was considered to indicate a statistically significant difference.

Evaluation of tumor-infiltrating immune cells

As described previously (Liu et al., 2017), the CIBERSORT method was used to quantify the proportions of 22 TIICs both in GC samples and normal samples using the LM22 signature and 1,000 permutations. Cases with CIBERSORT P < 0.05, which reflected that the deconvolution results were accurate, would be selected for further analysis. In the present study, a total of 222 samples (15 normal samples and 207 GC samples in TCGA dataset; 52 GC samples in the ICGC dataset) were filtered out. Group comparisons of 22 TIICs proportions were performed by t-test between normal and GC samples, and one-way ANOVA among different TNM stages using SPSS 20.0 (IBM Corp.). P < 0.05 was considered to indicate a statistically significant difference.

Construction and evaluation of a nomogram

LASSO Cox regression analysis was used to select TIICs that were highly associated with the overall survival of GC patients among the 22 TIICs. Subsequently, candidate TIICs were subjected to multivariate Cox regression analysis to further screen out the best survival-related candidate TIICs and develop an immunoscore formula. On the basis of immunoscore, patients were divided into low- and high-immunoscore groups. Ultimately, a prognostic nomogram integrating immunoscore and clinical information was developed for predicting 3-year survival probabilities of GC patients. Meanwhile, a calibration plot with bootstrapping set to 1,000 resamples and receiver operating characteristic curve (ROC) analysis was performed to assess the predictive capacity of the prognostic nomogram by calculating the area under the curve (AUC). Similarly, ROC curve analyses of TNM stage alone, the ICGC set and the whole set were also performed to validate the prognostic nomogram.

Expression analysis of differentially expressed genes (DEGs) between low- and high immunoscore groups

The statistical software program R (version 3.5.2; R Core Team, 2018) and the Bioconductor package edgeR (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) were used to identify DEGs between the low- and high-immunoscore groups, with the criteria of —log(fold-change)—>1.5 and false discovery rate <0.05.

Functional enrichment analysis and protein-protein interaction (PPI) network construction

All DEGs between the low- and high-immunoscore groups were utilized for Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analysis with adjusted P < 0.01 as the threshold. Meanwhile, all DEGs were inputted into STRING (https://string-db.org) (Szklarczyk et al., 2019) to predict protein-protein interactions, with a confidence >0.9 as the cut-off criterion. Data of the PPI network were processed by Cytoscape and genes in significant modules (MCODE score >10 and number of nodes >10) were extracted from the PPI network as the most important targets of differing immunoscores and were used for further analysis.

Survival analysis

Kaplan–Meier analysis and a log-rank test were performed to construct a survival curve and assess the survival difference between the low- and high-immunoscore groups, and analyze the association of the proportions of 22 TIICs with overall survival. In addition, univariate and multivariate Cox regression models were used to determine independent prognostic factors. P < 0.05 was set as cut-off value.

Expression profile of immunomodulators

Immune checkpoint inhibitors have revolutionized cancer therapy and have been approved for various cancer treatments. In the present study, several key immunomodulators (LAG3, TIM3, NKG2A, VISTA, CTLA-4, IFNG, IL2, IL6, ICOS, ICAM1, TIGIT, PD-1 and PD-L1) were quantified both in normal samples and GC samples. The differences in expression of the immunomodulators between normal and GC samples, as well as low- and high-immunoscore groups, were compared by t-test.

Patient characteristics

After cases with a follow-up time <30 days and CIBERSORT P > 0.05 were removed, 274 samples were enrolled in the present study, including 15 normal gastric samples and 207 GC samples in TCGA set and 52 GC samples in the ICGC set. All 259 patients were diagnosed pathologically with GC. The detailed demographic and baseline characteristics of the 259 GC patients are presented in Table 1. There were no statistically significant differences in the clinical information, including age, sex, ethnicity, tumor stage, immunoscore and survival status, between the two datasets.

Evaluation of tumor-infiltrating immune cells

With the CIBERSORT method, the present study quantified the proportions of 22 TIICs both in GC samples and normal samples, including naïve B cells, memory B cells, plasma cells, naïve CD4 T cells, resting memory CD4 T cells, activated memory CD4 T cells, CD8 T cells, gamma delta T cells, T follicular helper cells, regulatory T cells, resting dendritic cells, activated dendritic cells, macrophages (M0, M1 and M2), resting NK cells, activated NK cells, resting mast cells, activated mast cells, monocytes, neutrophils and eosinophils (Fig. 1A). In normal gastric samples, the proportion of plasma cells (29.3%) was the most abundant among 22 TIICs, followed by resting memory CD4 T cells (17.5%) and CD8 T cells (13.8%). However, in GC samples, the highest proportion among the 22 TIICs was macrophages (25.3%), followed by resting memory CD4 T cells (16.2%) and CD8 T cells (9.7%). As presented in Fig. 1B, the percentage of plasma cells (P < 0.001) and monocytes (P = 0.001) in GC samples was significantly lower than that in normal gastric samples, whereas, the fraction of activated memory CD4 T cells (P = 0.002) and macrophages (M0, P < 0.001; M1, P < 0.001; M2, P < 0.001) was opposite.

Construction and evaluation of a prognostic nomogram

Using the LASSO Cox regression model (Fig. 2A), TIICs that were highly correlated with the OS of GC patients were determined and subjected to multivariate Cox regression analysis to identify the best survival-related TIICs. Ultimately, seven TIICs (CD8 T cells, activated memory CD4 T cells, gamma delta T cells, monocytes, macrophages M2, neutrophils and eosinophils) were filtered out and used to develop an immunoscore formula (Fig. 2B). Based on the immunoscore, GC patients were divided into low- and high-immunoscore groups. The distribution of immunoscore, survival status and expression profile of the seven TIICs of each patient are presented in Figs. 2C–2E. No significant difference was identified between different TNM stages (Fig. 2F).

By integrating clinical information and the immunoscore, a prognostic nomogram for predicting 3-year survival probabilities of GC patients was built (Fig. 3). The AUC of the prognostic nomogram and TNM stage alone in TCGA set was 0.772 and 0.591, respectively (Figs. 4B and 4E). In addition, the AUC of the ICGC set and whole set was 0.730 and 0.748, respectively (Figs. 4C and 4D). The calibration plot of the prognostic nomogram for 3-year survival probability is shown in Fig. 4A.

Functional enrichment analysis and PPI network

A total of 342 DEGs (257 upregulated and 85 downregulated genes) were determined between the low- and high-immunoscore groups (Fig. S1). The results of KEGG analysis revealed that 342 DEGs were predominantly involved in 13 pathways, among which, a few pathways were highly associated with the human immune system, such as ‘Cell adhesion molecules (CAMs)’, ‘Chemokine signaling pathway’ and ‘Cytokine-cytokine receptor interaction’ (Fig. 5A). In GO analysis (Fig. 5B; Table S1), 147 biological processes were related to 342 DEGs, including ‘T cell migration’ and ‘regulation of T cell migration’.

A PPI network of DEGs with confidence >0.9 was gathered from the online STRING database and then processed by Cytoscape. There were 86 nodes (60 upregulated and 26 downregulated genes) and 262 edges in the PPI network (Fig. 5C). With MCODE score >10 and number of nodes >10 as the cut-off values, a MCODE containing 13 genes (CCL20, CXCL10, SAA1, CCL5, CXCR6, CXCL13, CXCL11, CXCL9, ADCY5, PYY, TAS1R3, NMU and CHRM2) was extracted from the PPI network and used for further analysis (Fig. 5D). KEGG analysis revealed that 13 genes in MCODE were mainly enriched in the following four pathways: ‘Chemokine signaling pathway’, ‘Cytokine-cytokine receptor interaction’, ‘Toll-like receptor signaling pathway’ and ‘TNF signaling pathway’ (Fig. S2).

Survival analysis

Kaplan–Meier analysis was performed to determine the association between the OS of GC patients and the immunoscore, as well as the proportions of 22 TIICs. The results of survival analysis indicted that the prognosis of GC patients in the high-immunoscore group was significantly worse compared with those in the low-immunoscore group (P < 0.001; Fig. 6A). In addition, a low proportion of activated memory CD4 T cells (P = 0.002; Fig. 6B) and CD8 T cells (P = 0.034; Fig. 6C), along with a high proportion of M2 macrophages (P = 0.019; Fig. 6D) was negatively associated with a favorable outcome for GC patients. All results of the univariate Cox regression analysis in TCGA [hazard ratio (HR), 2.406; 95% confidence interval (CI) [1.521–3.807]; P < 0.001], the ICGC (HR, 6.671; 95% CI [2.017–22.063]; P = 0.002) and whole sets (HR, 2.858; 95% CI [1.870–4.369]; P < 0.001) suggested that the association of immunoscore with OS was significant (Figs. 7A–7C). Furthermore, the results of multivariate Cox regression analyses in three sets indicated that the immunoscore was an independent prognostic factor (TCGA set: HR, 2.515; 95% CI [1.585–3.992]; P < 0.001; ICGC set: HR, 7.719; 95% CI [2.249–26.493]; P = 0.001; whole set: HR, 3.033; 95% CI [1.977–4.656]; P < 0.001; Figs. 7D–7F).

Expression profile of immunomodulators

As presented in Fig. 8A, TIM3, CTLA4, INF-γ, IL6, ICAM1, TIGIT and PD-L1 were significantly increased in GC samples compared with normal samples, whereas VISTA was significantly decreased. In addition, when compared with the low-immunoscore group, LAG3, TIM3, NKG2A, CTLA4, INF-γ, IL2, ICOS, ICAM1, TIGIT, PD-1 and PD-L1 were all significantly downregulated in the high- immunoscore group (Fig. 8B).