Construction and Validation of an Immune-Related Risk Score Model for Survival Prediction in Glioblastoma

Background: As one of the most important brain tumors, glioblastoma (GBM) has a poor prognosis, especially in adults. Immune-related genes (IRGs) and immune cell infiltration are responsible for the pathogenesis of GBM. This study aimed to identify new tumor markers to predict the prognosis of patients with GBM.

Methods: The Cancer Genome Atlas (TCGA) database and ImmPort database were used for model construction. The Wilcoxon rank-sum test was applied to identify the differentially expressed IRGs (DEIRGs) between the GBM and normal samples. Univariate Cox regression analysis and Kaplan–Meier analysis was performed to investigate the relationship between each DEIRG and overall survival. Next, multivariate Cox regression analysis was exploited to further explore the prognostic potential of DEIRGs. A risk-score model was constructed based on the above results. The area under the curve (AUC) values were calculated to assess the effect of the model prediction. Furthermore, the Chinese Glioma Genome Atlas (CGGA) dataset was used for model validation. STRING database and functional enrichment analysis were used for exploring the gene interactions and the underlying functions and pathways. The CIBERSORT algorithm was used for correlation analysis of the marker genes and the tumor-infiltrating immune cells.

Results: There were 198 DEIRGs in GBM, including 153 upregulated genes and 45 downregulated genes. Seven marker genes (LYNX1, PRELID1P4, MMP9, TCF12, RGS14, RUNX1, and CCR2) were filtered out by sequential screening for DEIRGs. The regression coefficients (0.0410, 1.335, 0.005, −0.021, 0.123, 0.142, and −0.329) and expression data of the marker genes were used to construct the model. The AUC values for 1, 2, and 3 years were 0.744, 0.737, and 0.749 in the TCGA–GBM cohort and 0.612, 0.602, and 0.594 in the CGGA-GBM cohort, respectively, which indicated a high predictive power. The results of enrichment analysis revealed that these genes were enriched in the activation of T cell and cytokine receptor interaction pathways. The interaction network map demonstrated a close relationship between the marker genes MMP9 and CCR2. Infiltration analysis of the immune cells showed that dendritic cells (DCs) could identify GBM, while LYNX1, RUNX1, and CCR2 were significantly positively correlated with DCs expression.

Conclusion: This study analyzed the expression of IRGs in GBM and identified seven marker genes for the construction of an immune-related risk score model. These marker genes were found to be associated with DCs and were enriched in similar immune response pathways. These findings are likely to provide new insights for the immunotherapy of patients with GBM.

Introduction

Glioblastoma (GBM) is characterized by high mortality and accounts for 48.6% of primary malignant brain tumors (1). The incidence of GBM increases with age, with the highest incidence seen among people aged 75–84 years. The patients usually have an unfavorable prognosis, with a median survival time of 8 months [95% confidence interval (CI): 8–9] (1). Therefore, survival prediction is important and urgent for the treatment of patients with GBM.

In recent years, rapid developments in the field of bioinformatics have aided the exploration of molecular characteristics of cancer. Many new molecular markers and molecular characterization systems of GBM have emerged, which provide guidance for understanding the mechanism of progression and promote diagnosis and treatment (2). Studies related to mutations of isocitrate dehydrogenase and platelet-derived growth factor receptor (3), as well as promoter methylation of methylguanine methyltransferase show that considerable progress has been made in the genetic research of GBM biomarkers (4–6). However, favorable treatment options and effective prognostic markers remain insufficient. Therefore, clinicians and researchers need to identify novel GBM biomarkers to improve therapeutic accuracy.

Of late, advances in immunotherapy have proved that the immune microenvironment plays an important role in tumor biology and established the promising prospects of immunotherapy in the future (7, 8). The GBM microenvironment is usually immunosuppressive and consists of infiltrating immune cells, such as microglia, natural killer cells, and neutrophils. For example, the results of immunohistochemical analysis in literature (9) showed that more than 70% of the analyzed human glioma samples (n = 105) had significant neutrophil infiltration. And a higher degree of neutrophil infiltration was associated with higher malignant glioma grades (10, 11). However, the microenvironment lacks T cells and specific types of non-immune components, such as neurons, astrocytes, and tumor cells (7). Compared with most other types of tumors, the abundance of infiltrating T cells in glioblastoma is lower (8, 12). In addition, neurons, astrocytes and tumor cells are not the main types of infiltrating cells in GBM microenvironment, and they may play a role in GBM by direct synapses, pruning synapses and promoting synaptic formation (13, 14). Furthermore, which immune genes and immune cells are associated with the prognosis of GBM is yet to be elucidated. Thus, it is necessary to enhance our comprehension of the potential immunopathological mechanism of GBM progression.

In this study, RNA-Seq data were acquired from The Cancer Genome Atlas (TCGA) database and were matched with the list of immune-related genes (IRGs) in the ImmPort database. Next, the differentially expressed IRGs (DEIRGs) were identified, and a prognostic model was constructed via survival analysis and Cox proportional hazard regression analysis. The Chinese Glioma Genome Atlas (CGGA) data were used to validate the predictive value of the model. In addition, the biological functions and action pathways of seven marker genes in GBM were examined, and their correlation with 22 kinds of immune cells was analyzed. The results are expected to provide new insights to understand the pathogenesis of GBM and establish the strategies for immunotherapy.

Methods

Data Acquisition in GBM

The gene expression data and individual clinical information of GBM patients were derived from the TCGA and CGGA datasets. On excluding patient data that lacked information on the survival status and survival time, clinical information of 174 patients with TCGA and 657 patients with CGGA were obtained.

Acquisition of IRGs

IRGs were exported from the ImmPort database, including IMMUNE_RESPONSE and IMMUNE_SYSTEM_PROCESS gene lists. A total of 332 genes were assessed in all (Supplementary Table 1).

Variance Analysis

The data downloaded from the TCGA and CGGA databases based on log2 | FC | were first normalized. Then, the TCGA data were used to construct the model, and CGGA data were employed to test the predictive performance of the model in populations from different sources. Then, 56,763 expressed genes from 174 patients in TCGA were intersected with the list of immune-related genes, and 879 genes with the expression data were acquired. The “limma” package of R software was utilized to analyze the expression differences in the immune-related genes in TCGA-GBM patients, and the DEIRGs were screened out [the screening conditions were log2 | FC | > 2 and False Discovery Rate (FDR) < 0.05]. The “ggplot2” and “pheatmap” packages of the R software were used to visualize the final results.

Screening and Identification of Marker Genes

The gene expression data of DEIRGs were integrated with the clinical data of TCGA-GBM patients including the survival status and the survival time data. The “survival” package of R software was employed for univariate Cox proportional hazard regression analyses and Kaplan–Meier analysis. P < 0.05 in both cases was considered to indicate the statistical difference. Then, the “survival” package of R software was employed for multivariate Cox proportional hazard regression analysis to further explore the prognostic potential of the obtained genes. The stepwise forward regression method was employed for multivariate Cox regression analysis, and the Akaike information criterion (AIC) was applied to avoid over-fitting. The genes with the largest likelihood ratios and the lowest AIC values were selected as the marker genes.

Construction of Prognostic Model Based on Marker Genes

By combining the regression coefficient results of multivariate Cox regression analysis with the expression data of each marker gene, we established a risk score predictive model using the following calculation formula:

$\begin{array}{l}RiskScore=({\beta }_1\times gen{e}_{1}expression)+({\beta }_2\times gen{e}_{2}expression)\\ +\cdots +({\beta }_n\times gen{e}_{n}expression)\end{array}$

Where β corresponds to the regression coefficient.

Evaluation and Validation of the Prognostic Model

According to the abovementioned model, we obtained the risk score of each patient and the grouping results of GBM patients. The survival curves were created using the Kaplan–Meier (KM) method, and the logarithmic rank test was employed to evaluate the different survival rates between the two groups. P < 0.05 was considered to indicate statistical significance. The “survivalROC” package of R software was used to draw the time-dependent receiver operating characteristic (ROC) curves, while the area under the curve (AUC) was applied to assess the specificity and sensitivity of the model. One and two years were defined as the time nodes. Survival analysis and AUC values were performed on CGGA data to validate the prognostic prediction performance of the prognostic model in Chinese GBM patients.

Functional Enrichment Analyses

Gene Ontology (GO) analysis was applied to identify the common functions of DEIRGs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was employed to identify the significant pathways for gene enrichment. The “clusterProfiler” package of R software was employed to display the enrichment results, and the significance threshold was set to 0.05. The STRING website was used to visualize the interactions between genes.

Correlation Analysis of Immune Cell Infiltration

The CIBERSORT online site was referred for analyzing the differences in the infiltrate abundance of 22 immune cells between tumor samples and normal samples (15). Data from the Tumor IMmune Estimation Resource (TIMER) 2.0 database were used to analyze the relationship between marker genes and 22 types of immune cells in GBM patients (16). Data from the TIMER database were used to demonstrate the differential functions of individual immune cells between the tumor and normal samples (8, 17).

Statistical Analyses

The ROC curves were employed to determine the diagnostic differentiation of the DEIRGs. KM analysis was performed to estimate the survival time of the two subgroups, while a logarithmic rank test was employed to determine the difference in prognosis (18). Cox proportional hazard regression analysis was applied to distinguish significant differences (18). All statistical calculation processes were conducted in the statistical software environment R version 4.1.0 or Microsoft Excel 2019. All P-values were two-tailed. The specific P-values and the other statistical methods are described throughout the study.

Results

Screening and Identification of DEIRGs

Initially, a flow chart was drawn to depict the analysis process (Figure 1). The TCGA-GBM data included 169 tumor samples and 5 normal samples. Upon comparing the expression abundance of IRGs between the tumor and normal samples, 198 DEIRGs satisfying both log2 | FC | > 2 and FDR < 0.05 were screened. The heat map revealed that 153 genes were upregulated and 45 genes were downregulated in the tumor samples (Supplementary Table 2). The volcano map signified that the log2 | FC | and –log10 (FDR) of RNA-Seq in the TCGA dataset. The genes with FDR < 0.05 and log2 FC > 2 (< -2) were marked with red (blue) dots. Figure 2 presents the results of gene screening.

FIGURE 1

Figure 1. The flow chart of the study protocol.

FIGURE 2

Figure 2. Heatmap and volcano plot of 198 DEIRGs in TCGA-GBM patients. (A) Heatmap of 198 DEIRGs in TCGA-GBM patients. Blue and red indicate the genes of a lower and higher expression. (B) Volcano plot of 198 DEIRGs in TCGA-GBM patients. Blue and red, respectively, represent the significantly downregulated and upregulated genes. The names of genes in the dense regions are hidden.

Establishment of Immune-Related Risk Score

The patients enrolled in the study were divided into the high- and the low-expression groups according to the median single-gene expression level. The single-gene KM curves showed that the hazard ratios (HRs) of the two groups remained the same throughout the study. Next, the DEIRGs were included in the univariate Cox regression analysis and the KM analysis assessed by the logarithmic rank test. Only DEIRGs with P < 0.05 in both cases were considered as candidate overall survival-related (OS-related) genes. By combining these two analyses, 15 genes were obtained (Table 1; Figure 3A). Single-gene KM curves of these 15 genes are shown in Supplementary Figure 1. To improve the independent forecasting ability of the model, these genes were included in the subsequent multivariate Cox regression analysis. Seven genes (LYNX1, PRELID1P4, MMP9, TCF12, RGS14, RUNX1, and CCR2) were obtained and were integrated into the prognosis model as marker genes (Table 2; Figure 3B).

TABLE 1

Table 1. The results of the univariate Cox regression analysis.

FIGURE 3

Figure 3. Forest plots depicting the results of univariate and multivariate Cox analyses. (A) The results of 13 DEIRGs were screened out by the univariate Cox analysis (We removed two of the genes with a large confidence interval). (B) The results of 7 prognostic DEIRGs were screened out by following multivariate Cox analyses.

TABLE 2

Table 2. The results of the multivariate Cox regression analyses.

Among the seven marker genes used to calculate the risk score, LYNX1, MMP9, RGS14, and RUNX1 were considered as the prognostic pathogenic genes, with HR > 1, while PRELID1P4, TCF12, and CCR2 were considered as the prognostic protection genes, with HR <1. The regression coefficients and gene expression data of the marker genes were used to establish a prognostic model. The following formula provides a method for calculating the risk score:

$\begin{array}{l}RiskScore=0.041\times LYXN1-1.335\times PRELID1P4+0.005\\ \times MMP9-0.021\times TCF12+0.123\times RGS14\\ +0.142\times RUNX1-0.329\times CCR2\end{array}$

Based on the score formula, the risk score of each GBM was obtained. Next, the patients were divided into the high- and low-risk groups by taking the median value as the critical value.

Accuracy Assessment and Validation of the Prognostic Model

The time-dependent ROC curves were employed to assess the predictive effect of the prognostic model. The AUC values of the model were 0.744 at 1 year, 0.737 at 2 years, and 0.749 at 3 years in the TCGA cohort (Figure 4) and 0.612 at 1 year, 0.602 at 2 years, and 0.594 at 3 years in the CGGA cohort (Figure 5).

FIGURE 4

Figure 4. The time-dependent ROC curves in the TCGA training cohort. (A–C) ROC curves to evaluate the predictive accuracy of the model for OS at 1-, 2-, and 3-years in the TCGA training cohort.

FIGURE 5

Figure 5. The time-dependent ROC curves in the CGGA validation cohort. (A–C) ROC curves to validate the predictive accuracy of the model for OS at 1-, 2-, and 3-years in the CGGA validation cohort.

Furthermore, survival analysis was performed in the model building and validation cohorts. Figure 6A shows the grouping results of all the patients with GBM in TCGA. The heatmap and the scatter map alluded that certain differences existed in the expression of the marker genes and OS (Figures 6B,D). These results were validated in the CGGA cohort (Figures 7A,B,D). The KM curves demonstrated that the low-risk group had a significantly better prognosis than the high-risk group, with a P < 0.05 (Figure 6C). Similar results were obtained in the CGGA validation cohort too, with P < 0.05 (Figure 7C). These findings suggest that the model has a good distinguishing ability and can be used as a prognostic biomarker in Chinese patients with GBM. In the training cohort and the validation cohort, the model in the first 3 years was better in distinguishing between high-risk and low-risk patients. In addition, we analyzed the effect of risk score on 15-year prognosis in the validation cohort, which shows that the model has a certain popularizing significance. However, in GBM cases, the survival time of more than 10 years is <1%, and long-term survivors are a very rare group (19). Shorter progression-free interval, larger diagnostic age, and surgical resection may be associated with a sudden decline in the 10-year survival probability of GBM patients (19).

FIGURE 6

Figure 6. Survival analyses of the prognostic model in the TCGA-GBM cohort. (A) Risk score scatter map. (B) Heat map of differential expression of 7 marker genes. (C) The KM curves of OS in the high- and low-risk groups. (D) Survival status scatter map.

FIGURE 7

Figure 7. The results of survival analysis of the prognostic model in the CGGA-GBM cohort. (A) Risk score scatter map. (B) Heat map of the differential expression level of 7 marker genes. (C) KM curves of OS in the 2 risk subgroups. (D) Survival status scatter map.

Enrichment Analysis

To infer the potential functions of DEIRGs in the carcinogenesis and development of GBM, the enrichment of GO terms and KEGG pathways of 198 DEIRGs were analyzed. GO enrichment analysis revealed 637 significantly enriched items. Supplementary Table 3 shows 26 GO items with the highest gene counts. T cell activation was found to be the most enriched in 584 items of biological processes, the outer surface of the plasma membrane was the most enriched in 6 items of cellular components, and receptor-ligand activity and signal receptor activator activity were the most enriched in 47 items of molecular function (Figure 8A). KEGG pathway enrichment analysis implied that the most significant pathway was the cytokine receptor interaction pathway (Figure 8B). Therefore, the DEIRGs mainly affect the internal functions of the immune cells, which may alter the immune response activity of the cells and affect the results of immunotherapy in the patients.

FIGURE 8

Figure 8. The results of functional enrichment analysis of 198 DEIRGs. (A) The results of GO analysis demonstrated the most-enriched biological functions. (B) The results of KEGG analysis show the most enriched biological processes. The X-axis represents the enrichment fraction, while the Y-axis represents the GO term and the KEGG pathway. The color indicates the P-value, while the circle size of the bubble chart represents the gene counts.

In addition, gene–protein interaction was explored on the STRING website (Supplementary Table 4). The website contains comprehensive known and predicted protein–gene interactions (20, 21). STRING was employed to explore the interaction between the marker genes and the realization mechanism of the regulatory functions. Figure 9 shows the interactions of the five marker genes.

FIGURE 9

Figure 9. The interaction network map of 5 marker genes. The network nodes represent the protein types, the edges represent protein–protein associations, and the line thicknesses represent the strength of the data supporting the interaction between the genes.

The interaction network map demonstrated that MMP9 and CCR2 were enriched in the central region of the network and exhibited strong interactions with several DEIRGs. In addition, RGS14 was linked to the central genes via interactions with ITGB2 and the RGS gene family. RUNX1 and TCF12 were connected to the central genes chiefly via interactions with the associated genes IL2RA and IL10.

Analysis of Tumor Immune Infiltration

Clinical studies on immunotherapy have confirmed that tumor-infiltrating immune cells can serve as a reference for the prognosis and immunotherapy of some solid tumors (22, 23). The tumor microenvironment plays a key role in tumor pathogenesis (24, 25). The CIBERSORT algorithm was applied to analyze the differences in the 22 immune cells between the tumor samples and the normal samples (Figure 10). The expression data of the 879 IRGs in the patients with GBM in TCGA were included, and the perturbation times were set to 1,000.

FIGURE 10

Figure 10. The violin plot shows the infiltration levels of 22 types of immune cells (blue and red represent normal and GBM samples, respectively).

The violin plot indicated that the M2 macrophages and the resting dendritic cells (DCs) had high invasion scores in the tumor samples (infiltration fraction > 0.2). This finding had implications for the important roles of these two types of cells in the tumor microenvironment. In the GBM tissues, tumor-associated macrophages mainly originated from the peripheral blood monocytes, and M2 polarization was induced by ARS2/MAGL signals from the tumor stem cells (26, 27). DCs can not only induce tumor antigen specific immune response by recognizing and presenting immature antigens, but also directly or indirectly participate in the regulation of immune system (28). In addition, there was a significant infiltration-level difference in the DCs between the normal and GBM tissues (P < 0.05). The accumulation of eosinophils has been found in a variety of central nervous system diseases, but the involvement of eosinophils in the immune response in GBM has only been preliminarily recognized in some case studies (29). Eosinophils may regulate the function of other immune cells by producing cytokines and chemokines and activate infiltrating immune cells in the tumor microenvironment (30). This may provide an explanation for GBM identification and immunotherapy.

Next, the TIMER 2.0 database was used to describe the correlation between the six marker genes and the DCs (Figures 11A–F). Furthermore, the TIMER database was used to generate KM curves for the DCs based on the TCGA-GBM gene expression data (Figure 11G).

FIGURE 11

Figure 11. The correlation between 6 marker genes and DC. (A) LYNX1. (B) MMP9. (C) TCF12. (D) RGS14. (E) RUNX1. (F) CCR2. (G) Kaplan–Meier analysis of the DC infiltrating levels from the TCGA database.

Figure 11 demonstrates the correlation coefficients ρ and P-values between the expression levels of six marker genes and the infiltration abundance of DCs. All ρ > 0 indicated that the six marker genes were positively correlated with DCs. Moreover, the ρ > 0.3 and P < 0.05 of LYNX1, RUNX1, and CCR2 suggested their significant positive correlation with DCs. These three genes may affect the tumor immune microenvironment and regulate the growth and invasion rate of the tumor cells by affecting the invasion of the DCs. In particular, Dendritic Cell-Specific Intercellular adhesion molecule-3-Grabbing Non-integrin (DC-SIGN) is highly expressed on immature DCs (31). And the activation of DC-SIGN promotes the transcription of MMP9 (32). Literature (33) shows that human monocyte-derived DCs express a large amount of RGS14, and the expression level of RGS protein family can be changed by toll-like receptor signals. Literature (34) found that TNF- α induces DC-SIGN expression in renal tubular epithelial cells and may be regulated through MTOR-RUNX1 pathway. In Helicobacter pylori infection model, the loss of CCR2 signal leads to the defect of differentiation and maturation of DCs, which affects the immune response function of DCs (35). Furthermore, the logarithmic rank P < 0.05 of DC was obtained. DCs can recognize and present foreign antigens and regulate the growth of tumors in the central nervous system by increasing the activity of natural killer cells and natural killer T cells (36). However, a small amount of circulating DCs in the body may inhibit endogenous tumor immunity (36). Individualized DC vaccines based on tumor-associated antigens have been developed on this basis and have been proved to be safe and effective in some GBM treatments (37). This finding signifies that the DCs are an important type of tumor-infiltrating immune cells equipped with the ability to identify GBM tumors.

Discussion

As one of the most invasive brain tumors, GBM is characterized by an unfavorable prognosis. Immunotherapy may prolong the survival time of the patients to some extent in the future. By analyzing the IRGs in the TCGA-GBM and ImmPort databases, 198 DEIRGs were identified. By employing several dimensionality reduction methods, such as difference analysis and Cox proportional hazard regression analysis, seven genes capable of predicting OS were screened, namely, LYNX1, PRELID1P4, MMP9, TCF12, RGS14, RUNX1, and CCR2. According to the Cox proportional hazard regression coefficients of the seven marker genes, a method for calculating the immune-related risk score was determined. On the basis of the sample coding rules of TCGA and CGGA databases (38, 39), patients with a risk score higher than the median value were considered to have a high risk for poor prognosis. The immune-related risk score can be used as a new prognostic predictor to explain the toxic role in immunobiology from cellular and molecular pathways (40). It may prevent the occurrence of immune-related adverse events by early identification of patients at risk of immunotherapy, so as to continue to maintain and give full play to the advantages of tumor immunotherapy (40). In addition, differential expression analysis of the marker genes revealed the significantly elevated expression of five genes (PRELID1P4, MMP9, TCF12, RUNX1, and CCR2) and the significantly lowered expression of two genes (LYNX1 and RGS14). Functional analysis showed that these genes were primarily involved in T cell activation and that they made a difference via cytokine receptor interaction pathways. In addition, the marker genes can be used as characteristic immune cell-specific genes, and they are also the expression products of some stromal cells in GBM microenvironment (8). This not only provides an explanation for the pattern of GBM-specific immune infiltration, but also opens a new window for the pathogenesis of immune infiltration in GBM.

The characteristics of the seven marker genes provide rich biological knowledge and potential therapeutic information about GBM. MMP9, one of the marker genes, is a oncogene of GBM (41). Studies have shown that the expression of MMP9 is upregulated in the GBM tissues and that it has the potential to induce the proliferation of GBM, thus worsening the prognosis of patients with GBM (41). In GBM patients, MMP9 is closely associated with hematopoietic progenitor cells, which may contribute to the development of specific therapies aimed at reducing HPC (42). Moreover, the association and interaction between heat shock protein 27 and MMP9 may contribute to the development of drugs that inhibit the infiltration and migration of GBM. This interaction also provides a new field of vision for the treatment of GBM (43). TCF12 has been proved to be an important target for miR-154 and can regulate the epithelial-mesenchymal transformation of GBM (44). P53/miR-154/TCF12 pathway may be involved in inhibiting the growth and invasion of GBM cells, so it can be used as a potential therapeutic target (44). Another marker gene, RGS14, is highly expressed in the caudate nucleus of the brain and thymus (45–47). Previous studies have shown that RGS14 targets the shuttle between centrosome and nucleoplasm and participates in the regulation of stress-induced cellular response (46). However, the role of RGS14 in the tumor tissues is yet to be elucidated. It was found that RGS14 had a low expression in GBM tissues, but it was potentially tumorigenic. More work is needed to discover the different roles of RGS14 in the GBM tissues and normal tissues. Another marker gene, RUNX1, is highly expressed in GBM, especially of the mesenchymal subtype (Mes). RUNX1 may promote the invasion of GBM cells by inducing TGF β signal transduction. The Mes GBM is widely seen in samples showing a high degree of necrosis. As an important transcriptional regulator of mesenchymal transformation, RUNX1 is also closely related to tissue necrosis (48, 49). This finding indicates that RUNX1 is a key signal node and that it might have a potential role in targeted therapy (50). Furthermore, CCL2 and its related receptor (CCR2) play an important role in brain tumors and are involved in regulating the migration of monocytes into the vascular endothelium. CCR2 inhibition can reduce tumor myeloid cells, suggesting that it may play a role in delaying the progression of gliomas (51–53).

The results of GO analysis showed that the most significant enrichment of biological process was T cell activation. This result suggests that immune-related genes may cause alterations in the morphology and maturity of T cells. Studies have revealed that PD-L1 present on the surface of GBM promotes the activation of PD-1 receptor in microglia, which inhibits the continuous proliferation of T cells and downregulates the cytotoxic activity of lymphocytes (54). Furthermore, regulation of leukocyte cell–cell adhesion was also highly enriched. For example, proteins regulated by IL-1 β may alter the frequency and speed of monocytes adhering to basal cells (55). The most significant enrichment of Cellular Component showed that the differential immune genes mainly performed biological functions on the extracellular side of the plasma membrane. The most significant enrichments of Molecular Function were the activities of receptor-ligand and signal receptor activator. For example, PD-1 ligand and PD-L1 receptor can negatively regulate T cell response and maintain homeostasis of the tumor microenvironment. Related studies have shown that PD-L1 can be used as a biomarker to assess the World Health Organization classification of GBM (54). For example, mutations or overexpression of epidermal growth factor receptors can activate the downstream signaling pathways, which may explain the pathogenesis of GBM (56, 57). On the other hand, the activity of the signal receptor activator may be correlated with the activation of the signaling pathway. CMTM6 is a key regulator of PD-L1, regulating T cell activity in vivo and in vitro (58). High levels of CMTM6 may inhibit the anti-tumor immune response of T cells in GBM, so the overexpression of CMTM6 is associated with poor prognosis and short overall survival of GBM patients (58). As for KEGG enrichment results, cytokine receptor interaction was the most important pathway. Moreover, other pathways, such as the interaction of viral proteins with cytokines and cytokine receptors and the JAK-STAT signaling pathways, were correlated with tumorigenesis (59, 60).

Furthermore, studies have shown that high immune infiltration in GBM predicts a poor outcome (61). Resting DC infiltration is a prominent feature of the GBM microenvironment in preclinical models and clinical samples. DC is an important immune cell in both innate and adaptive immune systems and has a complex array of functions and phenotypes in GBM microenvironment. Improvements in immunotherapy for DCs are in the experimental stage (62–64). In addition, the antitumor immune response of the DCs in the GBM tumor microenvironment may be realized by interaction with microglia, tumor cells, and T cells (62). Specifically, the combined application of CCR2 + HSCs and anti-PD-1 has been proven to significantly increase the median and long-term survival rates in the GBM model (65). The main mechanism of MMP9 in the brain is that it promotes enzyme activity and activates various cytokine- and chemokine-related immune/inflammatory reactions. Moreover, MM9 can promote leukocyte extravasation into the brain parenchyma and aid in destroying the blood–brain barrier. Its improper release may cause the development of brain tumors (66).

The TCGA-GBM-based risk score prognostic model was further validated in the Chinese samples, which indicated that our model has a certain degree of robustness. However, the results of this study are only based on bioinformatics and can only be used as a reference for preclinical research. In the future, more experiments are needed to prove it, such as flow cytometry and immunohistochemistry, to open up new ideas for the development of immunotherapy for GBM.

Conclusion

In this study, an immune-related risk score was established based on TCGA-GBM and ImmPort databases, and the predictive effect was validated using the CGGA-GBM database. Based on the findings, it could be suggested that the seven-gene risk score can independently evaluate the survival of patients with GBM. In addition, this score may partly reflect the immune microenvironment of tumor infiltration via the association between marker genes and DCs. The score is also closely related to the response rate of immunotherapy, thus serving as an immunotherapy reference for patients with GBM.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

WR and ZL contributed to research concepts and design, while WR and WJ contributed to data collection and statistical analysis. All authors contributed to the critical revision of important intellectual content of the manuscript.

Funding

This work was jointly supported by the National Natural Science Foundation of China (No. 81973560) and the Key Laboratory of TCM Encephalopathy of Zhejiang Province (Grant No. 2020E10012).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fneur.2022.832944/full#supplementary-material

References

1. Ostrom QT, Patil N, Cioffi G, Waite K, Kruchko C, Barnholtz-Sloan JS. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2013-2017. Neuro Oncol. (2020) 22 (12 Suppl. 2):iv1–iv96. doi: 10.1093/neuonc/noaa200

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Lee E, Yong RL, Paddison P, Zhu J. Comparison of glioblastoma (GBM) molecular classification methods. Semin Cancer Biol. (2018) 53:201–11. doi: 10.1016/j.semcancer.2018.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Yan H, Parsons DW, Jin G, McLendon R, Rasheed BA, Yuan W, et al. IDH1 and IDH2 mutations in gliomas. N Engl J Med. (2009) 360:765–73. doi: 10.1056/NEJMoa0808710

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Brat DJ, Verhaak RG, Aldape KD, Yung WK, Salama SR, Cooper LA, et al. Comprehensive, integrative genomic analysis of diffuse lower-grade gliomas. N Engl J Med. (2015) 372:2481–98. doi: 10.1056/NEJMoa1402121

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Network TC. Corrigendum: comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature. (2013) 494:506. doi: 10.1038/nature11903

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Reifenberger G, Wirsching HG, Knobbe-Thomsen CB, Weller M. Advances in the molecular genetics of gliomas - implications for classification and therapy. Nat Rev Clin Oncol. (2017) 14:434–52. doi: 10.1038/nrclinonc.2016.204

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Lim M, Xia Y, Bettegowda C, Weller M. Current state of immunotherapy for glioblastoma. Nat Rev Clin Oncol. (2018) 15:422–42. doi: 10.1038/s41571-018-0003-5

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, et al. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. (2016) 17:174. doi: 10.1186/s13059-016-1028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

9. De Leo A, Ugolini A, Veglia F. Myeloid cells in glioblastoma microenvironment. Cells. (2020) 10:18. doi: 10.3390/cells10010018

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Rahbar A, Cederarv M, Wolmer-Solberg N, Tammik C, Stragliotto G, Peredo I, et al. Enhanced neutrophil activity is associated with shorter time to tumor progression in glioblastoma patients. Oncoimmunology. (2016) 5:e1075693. doi: 10.1080/2162402X.2015.1075693

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Schernberg A, Nivet A, Dhermain F, Ammari S, Escande A, Pallud J, et al. Neutrophilia as a biomarker for overall survival in newly diagnosed high-grade glioma patients undergoing chemoradiation. Clin Transl Radiat Oncol. (2018) 10:47–52. doi: 10.1016/j.ctro.2018.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Oble DA, Loewe R, Yu P, Mihm MC Jr. Focus on TILs: prognostic significance of tumor infiltrating lymphocytes in human melanoma. Cancer Immun. (2009) 9:3.

PubMed Abstract | Google Scholar

13. Radin DP, Tsirka SE. Interactions between tumor cells, neurons, and microglia in the glioma microenvironment. Int J Mol Sci. (2020) 21:8476. doi: 10.3390/ijms21228476

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Pasqualini C, Kozaki T, Bruschi M, Nguyen THH, Minard-Colin V, Castel D, et al. Modeling the interaction between the microenvironment and tumor cells in brain tumors. Neuron. (2020) 108:1025–44. doi: 10.1016/j.neuron.2020.09.018

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. (2020) 48:W509–14. doi: 10.1093/nar/gkaa407

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Li T, Fan J, Wang B, Traugh N, Chen Q, Jun S, et al. TIMER: A web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. (2017) 77:e108–10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Abstracts from the 2021 Society for neuro-oncology's 26th annual scientific meeting education day. Neuro Oncology. (2021) 23 (Supplement_6):NP-NP. doi: 10.1093/neuonc/noab196

CrossRef Full Text | Google Scholar

19. Tykocki T, Eltayeb M. Ten-year survival in glioblastoma. A systematic review. J Clin Neurosci. (2018) 54:7–13. doi: 10.1016/j.jocn.2018.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. (2021) 49:D605–d12. doi: 10.1093/nar/gkab835

PubMed Abstract | CrossRef Full Text | Google Scholar

21. von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. (2003) 31:258–61. doi: 10.1093/nar/gkg034

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bremnes RM, Busund LT, Kilvær TL, Andersen S, Richardsen E, Paulsen EE, et al. The role of tumor-infiltrating lymphocytes in development, progression, and prognosis of non-small cell lung cancer. J Thorac Oncol. (2016) 11:789–800. doi: 10.1016/j.jtho.2016.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Badalamenti G, Fanale D, Incorvaia L, Barraco N, Listì A, Maragliano R, et al. Role of tumor-infiltrating lymphocytes in patients with solid tumors: can a drop dig a stone? Cell Immunol. (2019) 343:103753. doi: 10.1016/j.cellimm.2018.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Wu T, Dai Y. Tumor microenvironment and therapeutic response. Cancer Lett. (2017) 387:61–8. doi: 10.1016/j.canlet.2016.01.043

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Hinshaw DC, Shevde LA. The tumor microenvironment innately modulates cancer progression. Cancer Res. (2019) 79:4557–66. doi: 10.1158/0008-5472.CAN-18-3962

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yin J, Kim SS, Choi E, Oh YT, Lin W, Kim TH, et al. ARS2/MAGL signaling in glioblastoma stem cells promotes self-renewal and M2-like polarization of tumor-associated macrophages. Nat Commun. (2020) 11:2978. doi: 10.1038/s41467-020-16789-2

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zhou W, Ke SQ, Huang Z, Flavahan W, Fang X, Paul J, et al. Periostin secreted by glioblastoma stem cells recruits M2 tumour-associated macrophages and promotes malignant growth. Nat Cell Biol. (2015) 17:170–82. doi: 10.1038/ncb3090

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Yang L, Guo G, Niu XY, Liu J. Dendritic cell-based immunotherapy treatment for glioblastoma multiforme. Biomed Res Int. (2015) 2015:717530. doi: 10.1155/2015/717530

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Curran CS, Bertics PJ. Eosinophils in glioblastoma biology. J Neuroinflammation. (2012) 9:11. doi: 10.1186/1742-2094-9-11

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Kita H. Eosinophils: multifaceted biological properties and roles in health and disease. Immunol Rev. (2011) 242:161–77. doi: 10.1111/j.1600-065X.2011.01026.x

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Geijtenbeek TB, Torensma R, van Vliet SJ, van Duijnhoven GC, Adema GJ, van Kooyk Y, et al. Identification of DC-SIGN, a novel dendritic cell-specific ICAM-3 receptor that supports primary immune responses. Cell. (2000) 100:575–85. doi: 10.1016/S0092-8674(00)80693-5

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Yuan M, Zhang X, Zhang J, Wang K, Zhang Y, Shang W, et al. DC-SIGN-LEF1/TCF1-miR-185 feedback loop promotes colorectal cancer invasion and metastasis. Cell Death Differ. (2020) 27:379–95. doi: 10.1038/s41418-019-0361-2

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Shi GX, Harrison K, Han SB, Moratz C, Kehrl JH. Toll-like receptor signaling alters the expression of regulator of G protein signaling proteins in dendritic cells: implications for G protein-coupled receptor signaling. J Immunol. (2004) 172:5175–84. doi: 10.4049/jimmunol.172.9.5175

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhou S, Zhang L, Feng D, Luo M, Xie R, Yang K, et al. The mTOR-RUNX1 pathway regulates DC-SIGN expression in renal tubular epithelial cells. Biochem Biophys Res Commun. (2019) 519:620–5. doi: 10.1016/j.bbrc.2019.09.042

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Sun X, Zhang M, El-Zaatari M, Huffnagle GB, Kao JY. CCR2 mediates Helicobacter pylori-induced immune tolerance and contributes to mucosal homeostasis. Helicobacter. (2017) 22:12366. doi: 10.1111/hel.12366

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wheeler CJ, Black KL. DCVax-Brain and DC vaccines in the treatment of GBM. Expert Opin Investig Drugs. (2009) 18:509–19. doi: 10.1517/13543780902841951

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Wang QT, Nie Y, Sun SN, Lin T, Han RJ, Jiang J, et al. Tumor-associated antigen-based personalized dendritic cell vaccine in solid tumor patients. Cancer Immunol Immunother. (2020) 69:1375–87. doi: 10.1007/s00262-020-02496-w

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Blum A, Wang P, Zenklusen JC. SnapShot: TCGA-analyzed tumors. Cell. (2018) 173:530. doi: 10.1016/j.cell.2018.03.059

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhao Z, Zhang KN, Wang Q, Li G, Zeng F, Zhang Y, et al. Chinese glioma genome atlas (CGGA): a comprehensive resource with functional genomic data from Chinese glioma patients. Genomics Proteomics Bioinformatics. (2021) 19:1–12. doi: 10.1016/j.gpb.2020.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Dhodapkar KM. Autoimmune complications of cancer immunotherapy. Curr Opin Immunol. (2019) 61:54–9. doi: 10.1016/j.coi.2019.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Wang R, Zhang S, Chen X, Li N, Li J, Jia R, et al. EIF4A3-induced circular RNA MMP9 (circMMP9) acts as a sponge of miR-124 and promotes glioblastoma multiforme cell tumorigenesis. Mol Cancer. (2018) 17:166. doi: 10.1186/s12943-018-0911-0

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Huizer K, Sacchetti A, Dik WA, Mustafa DA, Kros JM. Circulating proangiogenic cells and proteins in patients with glioma and acute myocardial infarction: differences in neovascularization between neoplasia and tissue regeneration. J Oncol. (2019) 2019:3560830. doi: 10.1155/2019/3560830

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Rajesh Y, Banerjee A, Pal I, Biswas A, Das S, Dey KK, et al. Delineation of crosstalk between HSP27 and MMP-2/MMP-9: A synergistic therapeutic avenue for glioblastoma management. Biochim Biophys Acta Gen Subj. (2019) 1863:1196–209. doi: 10.1016/j.bbagen.2019.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhu G, Yang S, Wang R, Lei J, Ji P, Wang J, et al. P53/miR-154 pathway regulates the epithelial-mesenchymal transition in glioblastoma multiforme cells by targeting TCF12. Neuropsychiatr Dis Treat. (2021) 17:681–93. doi: 10.2147/NDT.S273578

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Yin W, Tang G, Zhou Q, Cao Y, Li H, Fu X, et al. Expression profile analysis identifies a novel five-gene signature to improve prognosis prediction of glioblastoma. Front Genet. (2019) 10:419. doi: 10.3389/fgene.2019.00419

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Cho H, Kim DU, Kehrl JH. RGS14 is a centrosomal and nuclear cytoplasmic shuttling protein that traffics to promyelocytic leukemia nuclear bodies following heat shock. J Biol Chem. (2005) 280:805–14. doi: 10.1074/jbc.M408163200

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Gerber KJ, Squires KE, Hepler JR. Roles for regulator of G protein signaling proteins in synaptic signaling and plasticity. Mol Pharmacol. (2016) 89:273–86. doi: 10.1124/mol.115.102210

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Zhao K, Cui X, Wang Q, Fang C, Tan Y, Wang Y, et al. RUNX1 contributes to the mesenchymal subtype of glioblastoma in a TGFβ pathway-dependent manner. Cell Death Dis. (2019) 10:877. doi: 10.1038/s41419-019-2108-x

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Wang Q, Cai J, Fang C, Yang C, Zhou J, Tan Y, et al. Mesenchymal glioblastoma constitutes a major ceRNA signature in the TGF-β pathway. Theranostics. (2018) 8:4733–49. doi: 10.7150/thno.26550

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Cooper LA, Gutman DA, Chisolm C, Appin C, Kong J, Rong Y, et al. The tumor microenvironment strongly impacts master transcriptional regulators and gene expression class of glioblastoma. Am J Pathol. (2012) 180:2108–19. doi: 10.1016/j.ajpath.2012.01.040

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Flores-Toro JA, Luo D, Gopinath A, Sarkisian MR, Campbell JJ, Charo IF, et al. CCR2 inhibition reduces tumor myeloid cells and unmasks a checkpoint inhibitor effect to slow progression of resistant murine gliomas. Proc Natl Acad Sci USA. (2020) 117:1129–38. doi: 10.1073/pnas.1910856117

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Chen Z, Feng X, Herting CJ, Garcia VA, Nie K, Pong WW, et al. Cellular and molecular identity of tumor-associated macrophages in glioblastoma. Cancer Res. (2017) 77:2266–78. doi: 10.1158/0008-5472.CAN-16-2310

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Groblewska M, Litman-Zawadzka A, Mroczko B. The role of selected chemokines and their receptors in the development of gliomas. Int J Mol Sci. (2020) 21. doi: 10.3390/ijms21103704

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Litak J, Mazurek M, Grochowski C, Kamieniak P, Roliński J. PD-L1/PD-1 axis in glioblastoma multiforme. Int J Mol Sci. (2019) 20:5347. doi: 10.3390/ijms20215347

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Shen CK, Huang BR, Yeh WL, Chen CW, Liu YS, Lai SW, et al. Regulatory effects of IL-1β in the interaction of GBM and tumor-associated monocyte through VCAM-1 and ICAM-1. Eur J Pharmacol. (2021) 905:174216. doi: 10.1016/j.ejphar.2021.174216

PubMed Abstract | CrossRef Full Text | Google Scholar

56. An Z, Aksoy O, Zheng T, Fan QW, Weiss WA. Epidermal growth factor receptor and EGFRvIII in glioblastoma: signaling pathways and targeted therapies. Oncogene. (2018) 37:1561–75. doi: 10.1038/s41388-017-0045-7

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Abstracts Abstracts from the 23rd Annual scientific meeting education day of the society for neuro-oncology november 15 – 18 2018 New Orleans Louisiana. Neurooncology. (2018) 20 (suppl_6):NP-NP. doi: 10.1093/neuonc/noy148

CrossRef Full Text | Google Scholar

58. Guan X, Zhang C, Zhao J, Sun G, Song Q, Jia W. CMTM6 overexpression is associated with molecular and clinical characteristics of malignancy and predicts poor prognosis in gliomas. EBioMedicine. (2018) 35:233–43. doi: 10.1016/j.ebiom.2018.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Rahaman SO, Vogelbaum MA, Haque SJ. Aberrant Stat3 signaling by interleukin-4 in malignant glioma cells: involvement of IL-13Ralpha2. Cancer Res. (2005) 65:2956–63. doi: 10.1158/0008-5472.CAN-04-3592

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Park AK, Kim P, Ballester LY, Esquenazi Y, Zhao Z. Subtype-specific signaling pathways and genomic aberrations associated with prognosis of glioblastoma. Neuro Oncol. (2019) 21:59–70. doi: 10.1093/neuonc/noy120

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Iglesia MD, Parker JS, Hoadley KA, Serody JS, Perou CM, Vincent BG. Genomic analysis of immune cell infiltrates across 11 tumor types. J Natl Cancer Inst. (2016) 108:djw144. doi: 10.1093/jnci/djw144

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Srivastava S, Jackson C, Kim T, Choi J, Lim M. A Characterization of dendritic cells and their role in immunotherapy in glioblastoma: from preclinical studies to clinical trials. Cancers. (2019) 11:537. doi: 10.3390/cancers11040537

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Finocchiaro G, Pellegatta S. Immunotherapy with dendritic cells loaded with glioblastoma stem cells: from preclinical to clinical studies. Cancer Immunol Immunother. (2016) 65:101–9. doi: 10.1007/s00262-015-1754-9

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Fecci PE, Sampson JH. The current state of immunotherapy for gliomas: an eye toward the future. J Neurosurg. (2019) 131:657–66. doi: 10.3171/2019.5.JNS181762

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Flores CT, Wildes TJ, Drake JA, Moore GL, Dean BD, Abraham RS, et al. Lin(-)CCR2(+) hematopoietic stem and progenitor cells overcome resistance to PD-1 blockade. Nat Commun. (2018) 9:4313. doi: 10.1038/s41467-018-06182-5

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Vafadari B, Salamian A, Kaczmarek L. MMP-9 in translation: from molecule to brain physiology, pathology, and therapy. J Neurochem. (2016) 139 (Suppl. 2):91–114.doi: 10.1111/jnc.13415

PubMed Abstract | CrossRef Full Text | Google Scholar