1. Association with stromal score and immune score with survival analysis and clinical parameters
375 samples with GC and 32 normal samples were downloaded from the TCGA database. The stromal scores and immune scores of all samples were obtained by the ESTIMATE algorithm. And the score ranges were -1859.703 to 2072.280 and -1056.272 to 3124.198, respectively. What’s more, the ESTIMATE score was -2471.016 to 4868.812.
To analyze the correlation between the stromal/immune/estimate scores and the overall survival of the samples, all samples were divided into high and low score groups. The Kaplan–Meier survival curve showed that the low score group of the stromal score has a higher 5-year survival rate than the high score group of the stromal score (Figure 1A). It suggested that the content of stromal cells in the GC microenvironment was related to the survival of patients. Survival rates were not significantly different in the high and low score groups of the immune/estimate scores (Figure 1B, C).
By exploring the potential relationship between the stromal/immune/estimate scores and the and clinical parameters, we observed that grade, stage, and tumor infiltration depth (T) had significant differences in scores (Figure 1D, E). In terms of the median value of stage, stage II, III, IV in the stromal and estimate median scores were higher than stage I, with statistical significance. From the aspect of the median value of grade, grade 3 was higher than grade 2 in the stromal score, grade 3 was higher than grade 2 and grade1 in the immune score, and grade 3 was higher than grade 1 in the estimate score, all with statistical significance. For the median value of tumor infiltration depth (T), T2, T3, and T4 were higher than T1 in the stromal/immune/estimate scores, with statistical significance.
Taken together, the correlation analyses of stromal and immune scores with clinical parameters of GC revealed that both high stromal and immune scores were related to poor tumor differentiation and severe local invasion and that neither was significantly related to lymph nodes and distant metastasis, indicative of the greater influence of TME on primary tumor cells than metastasis.
2. The characteristics of the tumor microenvironment and immunophenotypes
375 samples were analyzed using 29 immune gene sets according to the ssGSEA algorithm to evaluate the immune characteristics of GC patients. Based on the hierarchical clustering algorithm, the samples were clustered into three categories: 32 samples were clustered into cluster 1, 213 samples were clustered into cluster 2, and 128 samples were clustered into Cluster 3 (Figure 2A). As shown in Figure 2B, the immune characteristics of Cluster 1, Cluster 2 and Cluster 3 were compared,the immune characteristics of the three groups from high level to low level were: cluster1 characteristics (Immunity High group)> cluster 2 characteristics (Immunity Medium group)> cluster 3 (Immunity Low group).
Next, according to the scores of the tumor microenvironment (TME) of each sample, the TME characteristics in the Immunity High/ Medium /Low groups were analyzed (Figure 2C). For the Immunity High group, StromalScore, ImmuneScore, EstimateScore, and TumorPurity were 717.300±601.444, 2543.920±287.686, 3261.219±784.388, and 0.465±0.101, respectively. The Immunity Low group had the opposite trend: StromalScore, ImmuneScore, EstimateScore, and TumorPurity were-478.587±674.297, 251.226±426.450, -227.361±969.231, and 0.833±0.076, respectively. The Immunity Medium group was somewhere between Immunity High group and Immunity Low group: StromalScore, ImmuneScore, EstimateScore, and TumorPurity were 388.558±722.886, 1425.995±500.151, 1814.553±1072.298, 0.636±0.119, respectively.
The characteristics of the TME in the three groups were statistically significant differences (Figure 2E–G). In terms of tumor purity, the Immunity Low group had the highest level, while the Immunity High group had the lowest. Correspondingly, the Immunity High group had the highest infiltration degree of stromal cells and immune cells, and the lowest infiltration degree of stromal cells and immune cells in the Immunity Low group. All in all, with the decrease of tumor cell activity, the immune cell activity increased.
3. Identification and functional enrichment analysis of differentially expressed genes between high- and low-score groups
To disclose the potential relationship between the stromal/immune scores and the gene expression profile of the GC samples,we exploredthe differential analysis of transcriptome data of GC patients (Figure 3A, B). As can be seen from the Venn diagram, there were 640 identical up-regulated genes and 120 identical down-regulated genes (Figure 3C, D). The volcano plot revealed the top10 up-regulated genes and the top 10 down-regulatedgenes (Figure 3E).
We performed functional enrichment analysis on the obtained 760 differentially expressed genes(640 up-regulated genes and 120 down-regulated genes). Gene Ontology (GO)functions were mainly enriched in leukocyte proliferation,lymphocyte proliferation,mononuclear cell proliferation,regulation of lymphocyte activation,T cell activation,and lymphocyte differentiation (Figure 3F, G), while Kyoto Encyclopedia of Genes and Genomes(KEGG) pathways were mainly enriched in viral protein interaction with cytokine and cytokine receptor, staphylococcus aureus infection, hematopoietic cell lineage, chemokine signaling pathway, cytokine-cytokine receptor interaction, cell adhesion molecules (CAMs),and phagosome (Figure 3H, I). These results suggest that differentially expressed genes are associated with immune function.
4. Identification and verification of the prognostic-related hub gene
We constructed a PPI network of genes with prognostic values using the STRING network to analyze the interrelationship between genes with prognostic values. In this PPI network, there were 188 nodes and 372 edges, with red nodes representing up-regulated genes and purple nodes representing down-regulated genes (Figure 4A). The hub genes in the PPI network were identified by analyzing the number of adjacent nodes of each gene. Figure 4B showed the top 30 hub genes from the PPI network. To obtain the prognostic-related genes, we performed a univariate COX analysis of 760 differentially expressed genes (Figure 4C). The result showed that the 26 genes related to prognosis were high-risk genes (hazard ratio, HR>1), unfortunately, indicating poor prognosis.
The intersection of PPI hub genes and prognostic-related genes was taken to obtain the hub prognostic gene (Figure 4D). Analyzing the expression of the CXCR4 gene in the samples, we found that CXCR4 was significantly different between the normal and tumor groups (P<0.05), it was down-regulated in the normal group and up-regulated in the tumor group (Figure 4E). At the same time, we conducted a paired difference analysis on the samples of CXCR4 and found that the level of CXCR4 expression was significantly different in the paracancer group and the tumor group (P<0.05, Figure 4F). Finally, we further verified the CXCR4 expression in the Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) database, and the results were consistent with the above results (Figure 4G).
5. Analysis of CXCR4 gene
According to the median value of CXCR4 gene expression level, GC patients were divided into high expression group and low expression group. The study found that the 5-year survival rate of the low expression group was higher than that of the high expression group (Figure 5A). Clinical studies of the CXCR4 gene showed that the expression of CXCR4 in GC patients is significantly different in age, grade, stage, T, and N(Figure 5B-H).
Furthermore, CXCR4-related KEGG pathways were analyzed by GSEA software. GSEA analysis demonstrated that the CXCR4 gene was associated with several immune-associated signaling pathways, including intestinal immune network for IgAproduction (NES = 2.119, P < 0.001), leukocyte transendothelial migration(NES = 2.188, P < 0.001), natural killer cell mediated cytotoxicity(NES = 2.120, P < 0.001), and T cell receptor signaling pathway (NES = 2.113, P < 0.001) (Figure 5I, J).
6. Correlation Between CXCR4 and Tumor Immune Cell Infiltration
The CIBERSORT algorithm was employed to determine the degree of infiltration of immune cells in GC samples (Figure 6A). Moreover, different correlation patterns among the immune cells were found in TCGA cohorts (Figure 6B). AS a result, T cells CD4 memory activated was negatively correlated with B cells naïve, T cells regulatory (Tregs), and T cells CD4+ (CD4+ T) memory resting, but positively correlated with Macrophages M1 and T cells CD8+ (CD8+ T). Also, the contents of immune cells varied among different immune groups (Figure 6C).
The expression level of CXCR4 was positively correlated with the infiltration degree of immune cells (Figure 6D). By analyzing the differences of 22 immune cells in the high expression group and low expression group of CXCR4, we found that CXCR4 gene expression was significantly correlated with B cells naïve, B cells memory, T cells gamma delta, Monocytes, and Macrophages M0(Figure 6E). As shown in Figure F, B cells naïve,B cells memory,CD8+ T,T cells gamma delta, Mast cells activated,and Eosinophils infiltration were positively correlated with the high expression group of CXCR4, while Neutrophils,Macrophages M0,Mast cells activated,and Dendritic cells (DC) activated infiltration reached an inverse outcome.
7. The Prognostic Implication of CXCR4-related Immunomodulator in gastric Cancer
We identified 20 immunoinhibitors (ADORA2A, BTLA, CD160, CD244, CD274, CD96, CSFIR, CTLA4, HAVCR2, IDO1, IL10, KDR, LAG3, LGALS9, PDCD1, PDCD1LG2, PVRL2, TGFB1, TGFBR1, and TIGIT), 37 immunostimulators (CXCR4, ENTPDI, HHLA2, ICOS, ICOSLG, IL2RA, IL6, lL6R, KLRCI, KLRKI, LTA, PVR, TMEM173, TMlGD2, TNFRSF13B, TNFRSF13C, TNFRSF17, TNFRSF18, TNFRSF25, TNFRSF4, TNFRSF8, TNFRSF9, TNFSF13B, TNFSF14, TNFSF18, TNFSF4, and ULBP1) and 16 MHC molecules(MHC molecule:B2M, HLA-A, HLA-B, HLA-C, HLA-DMA, HLA-DMB, HLA-DOA, HLA-DOB, HLA-DPA1, HLA-DPB1, HLA-DQAI, HLA-DQA2, HLA-DQBI, HLA-DRA, HLA-DRB1, and HLA-E) significantly associated with CXCR4 in GC(Figure 7A). The volcano plot showed the 10 top genes that were tightly correlated to these immunomodulators (Figure 7B).
To investigate the prognostic values of CXCR4-associated immunomodulators in GC, we entered these variables into univariate and multivariate Cox regression analysis. This method led to an optimal 2-gene prognostic signature in GC. And we constructed a 2-gene prognostic model. The Cox regression analysis and the biological functions of the genes are presented in Table 1. The risk scores were calculated by adding up the product of expression value and coefficient of each gene. The Kaplan–Meier survival curve elucidated that patients with low-risk scores had a higher 5-year survival rate than those with high risk (log-rank test, P < 0.001) (Figure 7C). The area under the curve (AUC) value of the risk score was 0.598(Figure 7D). The distribution of risk scores, survival statuses, and the risk genes expression profiles for GC was visualized in Figure 7E.Figure 7F, G showed the correlation between the risk score and the survival rate in TCGA-STAD in univariate and multivariate COX regression models. Moreover, univariate Cox regression demonstrated that the risk score was an independent predictor of prognosis in GC after adjusting for age, gender, stage, T, M, and N (HR = 2.119, 95% CI = 1.196–3.755, P =0.010).
Table 1
Gene
|
Full name of gene
|
Functions
|
Coefficient
|
Univariate analysis
|
|
|
|
|
HR (95% CI)
|
p-value
|
CD86
|
CD86 molecule
|
coreceptor activity; protein binding; signaling receptor activity; virus receptor activity
|
0.0437
|
1.06(1.00-1.12)
|
0.044
|
TNFSF18
|
TNF superfamily member 18
|
cytokine activity; identical protein binding; protein binding; signaling receptor binding; tumor necrosis factor receptor superfamily binding; tumor necrosis factor receptor superfamily binding
|
0.1559
|
1.22(1.02-.47)
|
0.033
|
Lasso Cox analysis |
Finally, we constructed a prognostic nomogram in GC to anticipate the individuals' survival probability by weighing risk score, stage, grade,T, N, M, age, and gender (Fig. 7H). We also used time-dependent ROC to evaluate the predictive discrimination of the nomogram.(Fig. 7I). The AUC values of the 1-year, 3-year, and 5-year were 0.59,0.58 and 0.54, respectively.