Microarray dataset GSE62254, used as a training cohort for its complete clinicopathological and survival information, was downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/). GSE26253 dataset, which included 432 GC patients with complete survival information, was used as a testing cohort. The inclusion and exclusion criteria were exhibited in Figure 1 . The patients diagnosed as IV stage after surgery were filtered by the criteria since DFS was the main endpoint in this study.

Gene expression profiles of GSE62254 and GSE26253 datasets downloaded from GEO database were analyzed by ssGSEA, which classifies gene sets with common physiological regulation, chromosomal localization and biological functions (26). Gene markers of 28 immune cells [activated CD4⁺ T cells (CD4⁺ Ta), activated CD8⁺ T cells (CD8⁺ Ta), activated dendritic cells (aDC), CD56^bright natural killer cells (CD56⁺ NK), CD56^dim natural killer cells (CD56⁻ NK), activated B cells (Ba), central memory CD4⁺ T cells (CD4⁺ Tcm), central memory CD8⁺ T cells (CD8⁺ Tcm), effector memory CD4⁺ T cells (CD4⁺ Tem), effector memory CD8⁺ T cells (CD8⁺ Tem), eosinophils, gamma delta T cells (γδT), immature B cells (Bi), immature dendritic cells (DCi), mast cells, myeloid-derived suppressor cells (MDSC), memory B cells (Bm), monocytes, natural killer cells (NK), natural killer T cells (NK T), neutrophils, plasmacytoid dendritic cells (pDC), macrophages, regulatory T cells (Tregs), follicular helper T cells (Tfh), type-1 T helper cells (Th1), type-17 T helper cells (Th17), and type-2 T helper cells (Th2)] were obtained from a previous study (27). On the basis of the expression of these markers, the infiltration levels of immune cell types were quantified by ssGSEA in the R package GSVA (28).

To select the most useful prognostic immune cells, the “glmnet” package in R was utilized to perform the COX regression analysis with LASSO algorithm (29). Eight immune cells were selected from candidate cells by LASSO algorithm, where the data were subsampled and the tuning parameters were determined according to the expected generalization error estimated from 10-fold cross-validation. The optimal cut-off values were evaluated based on the association between DFS and cell fraction in the training cohort using the “survminer” package, dividing the patients into low-IS and high-IS groups.

Following the multivariate Cox regression analysis for the selection of independent prognostic factors, IS and other clinical pathological characteristics were used to generate the nomograms and calibration plots by “rms” package in R. In this model, each factor was assigned a weight score based on the results of the multivariate Cox regression analysis. Calibration was used to assess the performance of the nomograms. Receiver operating characteristic (ROC) analysis was also performed to estimate the accuracy of the nomograms for survival prediction using the “survival ROC” package of R. In addition, C-index was calculated with “survival” package.

Group comparisons were performed for continuous and categorical variables using one-way ANOVA and the χ² test. Survival curves were constructed by the Kaplan–Meier method and compared by means of the log rank test. Hazard ratios for univariable analyses were calculated using a univariable Cox proportional hazards regression model. A multivariable Cox regression model with the enter method was used to determine independent prognostic factors. Correlations between the immunoscore and mRNA expression of genes were analyzed by means of Pearson’s correlation test. The sensitivity and specificity of the survival prediction based on the immunoscore were depicted by a time-dependent receiver operating characteristic (ROC) curve, with quantification of the area under the ROC curve using the “timeROC” package. All statistical tests were two-sided and P < 0.050 was considered statistically significant. Statistical analyses were conducted using R software and SPSS^® version 19.0 (IBM, Armonk, New York, USA).

As shown by the flowchart ( Figure 1 ), 273 patients from the GSE62254 dataset and 365 patients from the GSE26253 dataset with DFS and OS information were included after applying the exclusion criteria. The patients’ clinicopathological characteristics are detailed in Table 1 . Of the 273 patients in the training cohort, 186 (68.1%) were men and the median age was 64 (28–84) years. In addition, the median DFS was 35.7 months and OS was 53.3 months. Of the 365 patients in the testing cohort, 239 (65.5%) were men and the median age was 53 (23–74) years. The median DFS and OS were 60.4 and 69.6 months, respectively, in the testing cohort. Figure S1 shows the distribution of the immune cells in patients in the training cohort ( Figure S1A ) and the relationship between 28 TIICs ( Figure S1B ).

Based on the relationship between 28 TIICs and DFS, we then constructed an immune-cell model. The forest plot in Figure 2A calculated by univariate Cox analysis shows the association between each immune cell subset and DFS. In general, CD4⁺ Ta (HR = 0.21, P < 0.001), aDC (HR = 0.44, P = 0.047), Th17 (HR = 0.35, P = 0.013), and CD56^- NK (HR = 0.21, P < 0.001) are protective factors for DFS. On the contrary, mast cell (HR = 2.6, P = 0.019) and pDC (HR = 2.3, P = 0.038) are risk factors. Then through LASSO Cox regression model eight immune cells were selected to build the IS and the formula was as follows: IS = 0.20810997 × DC + 0.79762920 × Mast cell - 0.8729771 × CD4⁺ Ta - 0.37309914 × CD8⁺ Tem - 0.04459008×Th17 - 0.88370283 × CD56⁺ NK - 0.12005145 × Ba-0.39656416×Bm ( Figures 2B, C ). We calculated an IS for each patient based on their personalized levels of the eight cells. The predictive accuracy of the model in the training cohort was assessed by time-dependent ROC analysis at 1, 2, and 3 years where AUC values were 0.733, 0.779, and 0.784, respectively ( Figure 2D ).

Using the optimum cut-off value (IS = -0.65) obtained by the “survminer” package, patients in the training cohort were divided into high-IS and low-IS groups. Figure 3A shows the distribution of clinicopathological characteristics between high-IS and low-IS groups. The Kaplan-Meier curve in the training cohort revealed that patients in the low-IS group had longer DFS (P < 0.001) as well as OS (P < 0.001) compared with the high-IS group ( Figures 3B, C ). The 5-year DFS and OS for low-IS group were 28.2% and 60%, respectively and 14.1% and 22.5% respectively, for high-IS group. Similar results were observed in the testing cohort (P < 0.001 for both DFS and OS, Figures 3D, E and Figure S2 ). The ROC analysis of the testing cohort indicated that the model could predict the prognosis of GC patients accurately ( Figure 3F ). The multivariate Cox regression analysis showed that high-IS and TNM stage were independent prognostic indicators for both poor DFS and OS in either training cohort ( Table 2 , Table S1 ) or testing cohort ( Tables S2 and S3). Additionally, different distributions of immune cells between high-IS and low-IS groups are illustrated in Figure 3G . The infiltration of CD4⁺ Ta, CD8⁺ Ta, aDC, Th17, and CD56^- NK cells in the TME of low-IS patients was significantly greater than that in high-IS patients. These data demonstrated that the eight-immune cells model, IS, could precisely predict the DFS and OS of GC patients with surgery.

We also assessed the relationships between IS, the status of relapse and survival, and the distribution of the eight cells in the training and testing cohorts. Figure S3 shows that patients in the high-IS group had more recurrence and death events than among those in the low-IS group, which further verified the accuracy of the model.

So far, the TNM stage has been regarded as a powerful indicator to predict the outcomes of patients with cancers, however, variable prognosis of patients was observed in clinical practice due to the heterogeneity. Here, to examine the value of the model, we respectively performed stratified analysis of the patients in stages I, II, and III in both the training and testing cohorts. Consistently, DFS and OS were all much longer in low-IS patients with stages I, II, and III GC than high-IS patients in either the training or testing cohort ( Figures 4A–C , Figures S4A–C ). Cox regression and forest plots were used to demonstrate that high-IS was a risk factor for both DFS and OS in stages I, II, and III ( Figure 4D ). Moreover, Figure S5A demonstrates that IS showed a better prognostic accuracy for DFS in GC patients with surgery than the TNM stage in the training cohort. Meanwhile, compared with the TNM stage alone, combining the data with IS showed a better prognostic value for DFS and OS in two cohorts ( Figure S5 ). In conclusion, the prognostic value of IS is independent of TNM stage. For GC patients with surgery, the combination of IS and TNM stage for prediction of DFS and OS is superior to that of TNM alone.

For GC patients, Lauren subtype and the status of microsatellite are the two recognized indicators to predict the outcomes and the efficacy of treatments. Therefore, we identified the value of IS based on Lauren classification and the status of microsatellite by using stratified analysis. Figures S6A, B show that low-IS patients had longer DFS (all P < 0.001 for intestinal and diffuse subtypes) and OS (P = 0.00046 for intestinal subtype, P < 0.001 for diffuse subtype, respectively) than high-IS patients. However, this phenomenon was not observed in mixed subtype patients ( Figure S6C ). Forest plots revealed that high-IS was the risk factor for DFS and OS simultaneously in both intestinal and diffuse subtypes ( Figure S6D ). Similarly, for patients with MSI, DFS and OS of high-IS patients were significantly shorter than those of low-IS patients (all P < 0.001, Figure S7A ). Similar results were obtained in patients with microsatellite stability (MSS) on DFS (P < 0.001) and OS (P < 0.001) ( Figure S7B ). Meanwhile, high-IS was also a risk factor for DFS and OS regardless of the status of microsatellite stability ( Figure S7C ). Taken together, IS remains a statistically and clinically significant prognostic model when stratified by GC subtypes.

For patients with surgery, especially stages II and III, adjuvant chemotherapy is indispensable. In the current study, the patients’ clinicopathological characteristics distributing in the group accepting adjuvant chemotherapy were similar to that in the group without chemotherapy after surgery ( Figure S8 ). As shown in Figure 5A , adjuvant chemotherapy could improve DFS and OS simultaneously in the total cohort (P = 0.00047 for DFS). For patients in the low-IS group, adjuvant chemotherapy improved DFS (P = 0.0061), which was not observed in the high-IS group (P = 0.19; Figure 5A ). The results from subset analysis in stage I showed no difference in DFS between patients receiving adjuvant chemotherapy and not receiving adjuvant chemotherapy ( Figure 5B ). However, the subset analysis in stages II and III patients demonstrated that adjuvant chemotherapy significantly improved DFS in the low-IS group (P = 0.0041), but no significant effect was found in the high-IS group (P = 0.071; Figure 5C ). Stages II and III GC patients with low-IS would obtain longer DFS by receiving chemotherapy after surgery. To explain this phenomenon, we performed the relationship between IS value and the expression of immune checkpoint regulators and inflammatory mediators. Notably, IS value was shown to be negatively correlated with PD-L1 (P < 0.001, r = -0.19), CD40 (P < 0.001, r = -0.23), CD47 (P < 0.001, r = -0.42), CTLA4 (P < 0.001, r = -0.37), GZMB (P < 0.001, r = -0.52), Tim-3 (P < 0.001, r = -0.29), ICOS (P < 0.001, r = -0.30), IDO1 (P < 0.001, r = -0.53), LAG3 (P < 0.001, r = -0.55), and IFNG (P < 0.001, r = -0.56), whereas the interleukin family and TGF showed no statistical correlation ( Figure 5D ). Further, we compared the expression of these immune checkpoint regulators and inflammatory mediator between low-IS and high-IS groups. As Figure S9 shown, the expression of PD-L1 (P<0.001), CD40 (P=0.001), CD47 (P<0.001), CTLA-4 (P<0.0001), GZMB (P<0.0001), Tim-3 (P<0.0001), ICOS (P<0.0001), IDO1 (P<0.0001), LAG3 (P<0.0001), and IFNG (P<0.0001) in low-IS group were all higher than that in high-IS group. In summary, these results indicated that IS could identify candidate stages II and III patients with surgery who would benefit from adjuvant chemotherapy.

It is generally known that patients with stage II-III should receive adjuvant chemotherapy after surgery according to NCCN guidelines. Next, we investigated the effects of chemoregimen in the training dataset. Patients with stage II-III GC were divided into two groups, patients with Xeloda plus cisplatin (XP, 83 patients) and patients with chemoregimen based on fluorouracil (5-Fu, 44 patients), according to the adjuvant chemoregimen. As shown in Figure 6A , patients with low-IS had significantly longer DFS regardless of the adjuvant chemoregimen (P = 0.0021 and P = 0.00011 for the XP group and 5-Fu group, respectively). However, as patients were divided into low-IS and high-IS groups based on the IS classifier, we found that in high-IS group patients who received XP regimen after surgery had longer DFS compared with patients who received regimen 5-Fu ( Figure 6B , P = 0.042), while in low-IS group patients there was no difference between two chemoregimen ( Figure 6B , P = 0.9). Cox regression analysis notably demonstrated that patients with high-IS might benefit from XP rather than 5-Fu after surgery ( Figure 6C ).

To provide a quantitative method in clinical practice to predict the probability of 1-, 2-, and 3-year DFS in patients with stages II and III GC after surgery, we established two nomograms integrating clinicopathological factors and IS on the basis of multivariate Cox regression analysis ( Figure 7 ). The c-index values were 0.7 and 0.677 for the nomograms with adjuvant chemotherapy and without adjuvant chemotherapy after surgery respectively, indicating a satisfactory overlap with actual observations. The two nomograms based on IS could be used to predict the prognosis of patients with or without adjuvant chemotherapy in clinical practice.