Impact Factor 5.085 | CiteScore 5.4
More on impact ›

Original Research ARTICLE

Front. Immunol., 20 October 2020 |

Identification of Immune Cell Infiltration and Immune-Related Genes in the Tumor Microenvironment of Glioblastomas

Sicong Huang1†, Zijun Song2†, Tiesong Zhang1, Xuyan He3, Kaiyuan Huang1, Qihui Zhang4,5, Jian Shen1* and Jianwei Pan1*
  • 1Department of Neurosurgery, The First Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, China
  • 2The First Affiliated Hospital, Institute of Translational Medicine, Zhejiang University School of Medicine, Hangzhou, China
  • 3The First Affiliated Hospital, School of Public Health, Zhejiang University School of Medicine, Hangzhou, China
  • 4Department of Neurology, Dong Fang Hospital, Beijing University of Chinese Medicine, Beijing, China
  • 5Department of Clinical Neurosciences, University of Calgary, Calgary, AB, Canada

Glioblastoma (GBM) is one of the most prevalent malignant brain tumors with poor prognosis. Increasing evidence has revealed that infiltrating immune cells and other stromal components in the tumor microenvironment (TME) are associated with prognosis of GBM. The aim of the present study was to identify immune cells and immune-related genes extracted from TME in GBM. RNA-sequencing and clinical data of GBM were downloaded from The Cancer Genome Atlas (TCGA). Four survival-related immune cells were identified via Kaplan-Meier survival analysis and immune-related differentially expressed genes (DEGs) screened. Functional enrichment and protein-protein interaction (PPI) networks for the genes were constructed. In addition, we identified 24 hub genes and the expressions of 6 of the genes were significantly associated with prognosis of GBM. Finally, the genes were validated in single-cell sequencing studies of GBM, and the immune cells validated in an independent GBM cohort from the Chinese Glioma Genome Atlas (CGGA). Overall, 24 immune-related genes infiltrating the tumor microenvironment were identified in the present study, which could serve as novel biomarkers and immune therapeutic targets.


Glioblastoma (GBM) is the most common primary malignant brain tumor accounting for approximately 80% of all primary malignant brain tumors, and has a dismal prognosis and poor quality of life, with a median overall survival (OS) often < 1 year. Hereditary syndromes and ionizing radiation are the most common risk factors for GBM (1). The standard care of GBM is surgical resection followed by concomitant radiation therapy and chemotherapy with temozolomide (TMZ). Although multiple treatments have improved due to the development of gene therapy, immunotherapy, vaccine therapy, and others (2), therapeutic options for managing recurrence in GBM are limited. Immune checkpoint inhibitors (ICIs) such as anti-programmed cell death protein-1 (PD-1)/programmed death ligand-1 (PD-L1) and anti-cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) have been extensively studied for both primary and recurrent glioblastomas in medical research. However, most of the clinical studies for GBM based on ICIs and trials with vaccine therapies have been unsuccessful. The cause of the failure in clinical trials of GBM via immunotherapy is attributed to several factors, including a highly immunosuppressive environment and multiple mechanisms of therapeutic resistance. GBM induces local immune dysfunction and systemic immunosuppression, which causes more complex coupling relationships between GBM and the surrounding tumor microenvironment (TME). Studying the mechanisms of GBM immunosuppression enhances our understanding on development of immunotherapy strategies (3).

TME is one of the crucial factors of local immune dysfunction, which establishes a niche for cancer cells, multiple stromal cells (endothelial cells, immune cells, etc.) and extracellular components (extracellular matrix, cytokines, growth factors, etc.). TME plays a critical role in the establishment of specific conditions, thereby interfering with angiogenesis, cell death, oxidative stress, and immune escape (4). Increasing studies have revealed that TME is not only pivotal in tumor initiation, progression, and migration, but it also affects generation of therapeutic resistance and malignancy. Cellular composition of TME and accessibility of immune cells exhibit large variations among GBM subtypes and patients. Such factors contribute to immunosuppression of GBM, which in turn lead to immunotherapeutic treatment failure (5). Identification of actively involved types of immune genes and immune cells associated with the TME facilitates elucidation of the general mechanisms of GBM immunosuppression.

Therefore, the present study investigated survival-related immune cells in GBM and identified hub genes associated with immune cell infiltration. We acquired RNA-sequencing (RNA-seq) expression data and corresponding clinical data of 166 patients with GBM from The Cancer Genome Atlas (TCGA) database. A total of 22 types of infiltrating immune cells in the 166 patients were estimated using the method of estimating relative subsets of RNA transcripts (CIBERSORT) (6). Subsequently, four survival-related immune cells were identified from the survival analyses of 22 types of immune cells. Immune-related genes were ranked through differential gene expression analyses and 24 hub genes selected from the protein-protein interaction (PPI) network established using Cytoscape (7). Six hub genes associated with overall survival were identified. Finally, immune cells were validated in an independent GBM cohort from the Chinese Glioma Genome Atlas (CGGA), and hub genes verified in single-cell sequencing studies of GBM. All analyses were conducted using R software. The findings of the present study provide valuable information that will guide patient-specific clinical immunotherapeutic strategies, and further construction of prediction models for prognosis of GBM. Moreover, immune cells infiltrating in the tumor microenvironment could act as therapeutic targets for the clinical treatment of GBM.

Materials and Methods

Raw Data Collection

RNA-Seq expression profiles of immune cells and corresponding clinical data of 166 patients with GBM were downloaded from TCGA database. The file format of RNA-seq expression was FPKM. The expression profile of each sample included age, gender, expression subclass, and MGMT promoter status. RNA-Seq expression information of immune cells from CGGA were also downloaded for the validation. Data acquisition and analyses were performed using R software (8).The entire research data analysis process is presented in Figure 1.


Figure 1 Flow chart of the whole analysis process.

Identification of Survival-Related Tumor-Infiltrating Immune Cells

CIBERSORT is an analytical algorithm, which can characterize cell composition of complex tissues based on normalized gene expression profiles (9). We used CIBERSORT to estimate the ratio of 22 infiltrating immune cell types based on each GBM sample. Afterward, 57 samples with P ≤ 0.05 were selected and correlation analyses conducted to analyze contents of the 22 immune cells (10). Survival analyses of the filtered immune cells in the tumor microenvironment were performed by the Kaplan-Meier survival analysis, with a cut-off level set at the median value. The results were tested by log-rank test. All the analyses were conducted using R software.

Relationship Between Clinical Characteristics and Survival-Related Immune Cells

To determine the relationship between survival-related immune cells and clinical features such as age, gender, expression subclass, and MGMT promoter status, 57 samples were analyzed. An independent sample t-test was used to compare means of two groups, while one-way analysis of variance (ANOVA) test was used to compare the means of four groups.

Identification and Functional Enrichment Analysis of Immune-Related Genes

Immune related-genes were analyzed using survival-related cells that had been obtained previously. Data analysis was performed using the edgeR R package, and |logFC| ≥ 1.0 and P < 0.05 were set as the cut-offs to screen for immune-related genes. Subsequently, a Venn diagram was used to visualize genes displayed by the online tool ( (11). DAVID software ( was used to analyze immune-related genes in the Gene Ontology(GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (12). Results of GO analysis revealed the functions of immune-related genes in biology process, cell component, and molecular function (13). KEGG pathway analyses results revealed the role of development-related signaling pathways.

Construction of PPI Network, Selection and Analysis of Hub Genes

PPI networks of immune-related genes were predicted using the Search Tool for the Retrieval of Interacting Genes (STRING, (14). An interaction combined score of >0.4 was considered statistically significant. Cytoscape is an open-access software platform designed to analyze and visualize complex interaction networks (7). Molecular Complex Detection(MCODE) plugin of Cytoscape was used to cluster the networks based on topology to identify densely connected regions with MCODE score > 5, degree cut-off = 2, node score cut-off = 0.2, max depth = 100, and k-score = 2 (15). Hub genes were defined based on module connectivity (16).

Identification and Immune Infiltration of Survival-Related Hub Genes

Kaplan-Meier plots were used to identify immune-related genes in relation to the overall survival of patients. These results were analyzed by long-rank test. The correlation between 24 hub genes and 22 immune cells was determined using Person’s correlation analysis and CIBERSORT to reveal the relationship between hub genes and immune cells (17). Afterward, comprehensive correlation analysis between six selected survival-related hub genes and tumor-infiltrating immune cell signatures for GBM were performed using Tumor Immune Estimation Resource (TIMER 1.0, (18).

Distribution of Immune-Related Hub Genes in TME of GBM From Single-Cell Data

Data for the single cell GBM analysis was derived from the paper “An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma”, and the Seurat R package was used to reprocess the count matrix in which the dimensional reduction plot and cell type annotation were both retrieved from published meta data (19). The distribution of expressions of the hub genes was created using the Feature Plot function.


Data Source and Identification of Survival-Related Immune Cells

The workflow of the study is presented in Figure 1. Publicly available data for the 166 cases of GBM, including RNA-Seq (FPKM and counts format) and clinical data were downloaded from TCGA database. The abundance ratios of 22 immune cells in the 57 samples are presented in the Figure 2A, and the relationship between abundance ratios of the immune cells is presented in Figure 2B. Consequently, the correlations between abundance ratios of immune cells were analyzed using Kaplan-Meier survival analysis to elucidate the potential role of the abundance ratios of immune cells in overall survival. The four immune cells that were associated with survival are presented in Figures 2C–F. The results of survival analyses indicated that there was a significant negative correlation between M0 Macrophages, while monocytes, activated NK cells, and eosinophils predicted positive overall survival.


Figure 2 The abundance ratios of 22 immune cells and overall survival analysis. (A) The abundance ratios of immune cells in the 57 samples. The specific 22 immune cells corresponded to to one sample by different colors as shown in barplot. (B) The abundance ratios matrix of 22 immune cells. The value represents the correlation value, green represents the positive correlation while brown represents negative correlation. (CF) Overall survival analysis of four immune cells based on Kaplan Meier-plotter from the comparison of groups of high (yellow line) and low (blue line) genes expression. (p<0.05).

Clinical Data Correlated With Survival-Related Immune Cells

To determine the effect of immune cells on the clinical characteristics of GBM, relevant GBM clinical data were downloaded to investigate correlation with the abundance ratios of survival-related immune cells. The clinical characteristics included age, gender, expression subclass, and MGMT promoter status. The odds ratio of monocytes and eosinophils increased in neural and proneural types and was higher in males than in females (Figure 3).


Figure 3 Relationship between four survival-related immune cells and clinical features. (A–D) The relationship between four survival-related immune cells and age, gender, expression subclass, and MGMT status.

Screening of Immune-Related Genes

The immune-related genes were categorized into high- and low-expression groups in GBM to identify genes that associated with the four survival-related immune cells. Unique genes expression profiles of the four survival-related immune cells are presented by volcano plots in Figure 4. A total of 1,107 genes were identified in monocytes, 1,137 genes in macrophages M0, 1,742 genes in activated NK cells, and 1,336 genes in eosinophils (Figures 4A–D). In addition, 38 identical genes expressed in infiltration of the four immune cells are presented by Venn diagrams in Figure 4E.


Figure 4 Screening for immune-related genes. (A–D) The volcano plot of all quantified genes in the analysis of monocytes, macrophages M0, NK cells activated, and eosinophils. (E) Venn diagram indicates the overlap of differentially expressed genes across the four different immune cells.

Functional Enrichment Analysis of Immune-Related Genes

Functional enrichment analysis of immune-related genes was performed via DAVID website to reveal the potential functions of immune-related genes (Figure 5). GO term analysis revealed that immune-related genes were significantly enriched in the biological processes (BP) of nervous system development, cell adhesion, extracellular matrix organization, and chemical synaptic transmission (Figure 5A). Genes in the cellular components (CC) groups (Figure 5B) were primarily enriched in the plasma membrane, extracellular exosome, extracellular space, and extracellular region; the molecular functions (MF) were enriched in protein binding, calcium binding, structural constituent of cytoskeleton, and microtubule binding (Figure 5C). Moreover, the KEGG analysis revealed that immune-related genes were linked to cell adhesion molecules, cAMP signaling pathway, leukocyte transendothelial migration, protein digestion and absorption, and Toll-like receptor signaling pathway (Figure 5D). These results demonstrated that the genes were associated with the extracellular matrix of tumor microenvironment and cellular interaction.


Figure 5 Functional enrichment analysis of immune-related genes. (A) Biological process analysis. (B) Cellular components analysis. (C) Molecular function. (D) KEGG pathway analysis (p < 0.05).

Modular Analysis Based on PPI Network

Considering the limitation of the PPI networks regarding the number of genes, we screened all the differentially expressed genes but selected the genes only co-expressed in at least two immune cells. Overall, we identified 920 genes from 4,122 genes. These genes were imported into the online STRING tool to elucidate the interaction of immune-related genes. Finally, we got the PPI network with 357 genes which the combined-score was set to ≥ 0.4 (Figure 6A). We selected the most significant module for further functional enrichment analysis (Figure 6B).


Figure 6 Protein-Protein interaction network construction and modular analysis. (A) PPI network was constructed using a total of DEGs. (B) The most significant module was marked. The color of a node reflects the log(Fc) value of the gene expression, the size of a node suggests the numbers of interacting genes with others.

Identification of Survival-Related Hub Genes

A total of 24 hub genes with high connectivity in the modules were identified from the PPI network based on the cut-off criteria. We subsequently elevated the biological enrichment analysis of the 24 hub genes using the online tool ( (Table 1). Six of the hub genes were significantly correlated with survival (Figure 7). GRIA1, BST2, B2M, and TRIM21 were positively correlated with the overall survival. GRIA2 and MAP2 were correlated with poor prognosis. The relationship between 24 hub genes and 22 immune cells analyzed using Person’s correlation analysis is performed in Figure 8A. The remarkable relationship between infiltration levels of immune cell types and survival-related hub genes was validated in TIMER. The results indicated that infiltration levels of CD8+ T cells, neutrophils, and dendritic cells were significantly associated with GRIA1, GRIA2, and MAP2 (Figure 8B). Furthermore, BST2 and B2M were correlated with B cells, macrophages, and dendritic cells, and TRIM21 was associated with B cells and neutrophils.


Table 1 The function of hub genes.


Figure 7 Overall survival analysis of six hub genes. (A) B2M. (B) BST2. (C) GRIA1. (D) GRIA2. (E) MAP2. (F) TRIM21 (p < 0.05).


Figure 8 Immune infiltration of survival-related genes. (A) The correlation between expression proportion of hub genes and immune cells. Red suggests the positive correlation while the blue represents negative correlation. The size of point indicates P-value, and the color reflects the correlation. (B) The correlation analysis between survival-related genes and tumor infiltrating immune cells was performed. Scatter plots were generated with partial Spearman’s correlation and statistical significance.

Validation of the Correlation Between Immune Cell Infiltration and Survival-Related Hub Genes

The correlation between survival-related hub genes and immune cell infiltration in GBM was analyzed after determining the prognostic value of hub genes (Figure 9).


Figure 9 Validation of hub genes in single- cell sequencing in GBM. t-distributed neighbor embedding(tSNE) plot of all single cells. The color represents the expression of markers for Malignant cells (green), marcophages(magenta), oligodendrocytes (cyan), and T-cells (blue).

In addition, gene expression data of immune cells for 134 GBM samples were downloaded from CGGA database to investigate the significance of immune cells identified from TCGA database. The results we obtained from CCGA revealed that activated NK cells (Figure S1, p = 0.019) and monocytes (Figure S1, p = 0.023) were associated with positive prognosis, which are consistent with the data we have gotten previously from the TCGA database (Figure S1).

Validation of the Expression of Immune-Related Hub Genes by Single-Cell Sequencing

The cells were classified as malignant and non-malignant cell types by combining three approaches; high expression of markers classified as non-malignant cells such as macrophages, T cells, and oligodendrocytes. The distribution of hub genes expressions in the four cell clusters is displayed in Figure 9. A1F1, C3AR1, FCGR1A, MNDA, HMOX1, and TLR2 were only expressed in macrophages. B2M, CCT3, HSPA8, and TUBA1A were significantly expressed in all the four cell clusters. With reference to survival-related genes, BST2 was detected in macrophages, T-cells, and malignant cells. GRIA1 and GRIA2 were expressed in oligodendrocytes and malignant cells. MAP2 was only detected in malignant cells. However, TRIM21 was not detected in any of the cells types. Notably, microglia are the vital macrophages of the brain, and they act as the primary form of immune defense in the central nervous system. A specific microglial marker in humans, TMEM 119, was used to distinguish microglia from macrophages in the brain (Figure S2). We subsequently identified the expression of hub genes in microglia and found AIF1, B2M, BST2, C3AR1, CCND1, CCT3, FCGR1A, GNG7, HMOX1, HSPA8, MNDA, TLR2, and TUBA1A were significantly expressed (Figure S2).


The present study analyzed immune cells and immune-related genes in TME of GBM to establish a potential strategy for GBM immunotherapy. The study identified immune-related genes in TME, which significantly contributed to the survival of patients with GBM from TCGA database. Four survival-related immune cells were initially identified from GBM samples and the genes correlating to the levels of four immune cells analyzed. Furthermore, GO and KEGG enrichment analysis were conducted to investigate the biological functions of immune-related genes. Subsequently, all the immune-related genes were imported to construct a PPI network, and 24 hub genes obtained. Finally, the immune cell types in patients with GBM were validated using CGGA database, and hub genes validated in single-cell sequencing.

Four types of survival-related immune cells associated with GBM were identified from TCGA database, including M0 macrophage, monocytes, NK cells and eosinophils. Previous research has indicated that immune cells, especially tumor-associated macrophages (TAMs) in TME interact with tumor cells through direct contact or different signaling pathways. TAMs are crucial components of infiltrating immune cells, accounting for 30–40% of the cellular components in GBM (20). Immune cell populations in GBM are classified into two categories: microglia and bone marrow-derived monocytes. The BBB is damaged during tumor progression (21). With the accumulation of a family of monocyte chemoattractant family of proteins (MCPs), monocytes from the periphery infiltrate into the tumor across the BBB, and then differentiate into macrophages. Tumor-associated macrophages are often regarded as the facilitators of tumor proliferation due to their proangiogenic and immunosuppressive effects (21). M0 macrophages, which are referred to as ‘alternatively activated macrophages,’ can be polarized into M1 or M2 phenotypes by environmental signals (22). M1 macrophages can produce pro-inflammatory cytokines that are essential for host defense and exert tumoricidal effects in GBM (21). However, M2 macrophage phenotype is considerably involved in tumor cell proliferation and prediction of poor clinical prognosis in patients with GBM patients (23). M1 and M2 macrophages are plastic and heterogeneous immune cells, and the TME facilitates the regulation of functional polarization of TAMs (24). Currently, researchers have been working on promoting the reversal of TAMs from M2 to M1 based on their polarization (25, 26). Therefore, the results may indicate that the macrophages in TME of GBM could be used as potential therapeutic targets for GBM immunotherapy.

NK cells accounts for 2.11% of the total cellular components in GBM, which constitutes the lowest proportion of all immune cells infiltrating in GBM (27). NK cells have been reported to recognize target cells that are deficient in the surface expression of major histocompatibility complex (MHC) molecules, and can directly lyse tumor cells without prior activation (28). However, TME influences the immune function of NK cells and causes immune evasion. The upregulation of growth factor signaling pathways or the loss of cell cycle regulators promotes evasion of GBM from surveillance through resistance to NK-derived cytotoxicity (29). Moreover, GBM cells express high levels of MHC class I molecules and human leukocyte antigens (HLA)-A, HLA-B, and HLA-C ligands, which inhibit functions of NK cells via killer immunoglobulin-like receptors (KIRs) (30). Therefore, blocking KIRs could disrupt the tumor microenvironment and attenuate the activity of NK cells to kill GBM cells. Increasing the number of NK cells infiltrating the GMB microenvironment and modification of NK cells could be a potential treatment intervention for GBM (31, 32). Emerging evidence has demonstrated that the activation of eosinophils induces initiation, promotion and progression of GBM (33).

Previous advances have indicated the eosinophil-derived neurotoxin (EDN) and eosinophil cationic protein (ECP) play a critical role in preventing GBM initiation (34). During GBM promotion, eosinophils are activated by GBM mediators, which in turn lead to the production of tumor promoting growth factors (35). Nevertheless, the mechanisms of immune response in GBM remain indeterminate; therefore, further studies are required to investigate the mechanism involved.

More importantly, KEGG enrichment analysis indicated that these differential immune-related genes were enriched in the classical pathway, such as cell adhesion molecules (CAMs) and cAMP signaling pathway. CAMs are glycol-proteins expressed on the cell surface and play a critical role in multiple biologic processes during tumor development (36). It has been reported that CAMs mediate the process of immune responses in the tumor microenvironment, such as immune cell recruitment, immune cell activation, and formation of immunological synapse between immune cells and tumor cells (37). The cAMP signaling pathway, which acts as universal second messengers regulates pivotal physiological processes. The increases of intracellular cAMP inhibits innate immune functions (38). At the BP level, these differential immune-related genes were significantly enriched in cell adhesion, and extracellular matrix organization. In the CC groups, the differential immune-related genes were related to extracellular exosome, extracellular space, and extracellular region; the MF groups were enriched in protein binding, the structural constituent of cytoskeleton, and microtubule binding. Theses pathways are all related to the extracellular matrix components and cell’s cytoskeleton in the microenvironment. These above results further indicate the reliability of the immune differential genes and their relevance to the GBM tumor microenvironment.

Furthermore, we identified 24 hub genes, and 6 of these genes (GRIA2, GRIA1, BST2, MAP2, B2M, and TRIM21) have significant correlation with prognosis and were considered as predictive biomarkers that could provide valuable insights into new immunotherapy strategies. Previous studies have demonstrated that glioma cells can secrete excitotoxicity glutamate that mediates neuronal death in glioma microenvironment. Moreover, glutamate secretion promotes tumor expansion by inducing inflammatory response within the surrounding areas (39). Researchers have established that the expression of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptors (AMPARs) protects GBM cells from the glutamate-rich tumor microenvironment (40). AMPARs are complexes consisting of four subunits (GluR1, GluR2, GluR3, and GluR4). GRIA1 and GRIA2 are also referred to as GluR1and GluR2, respectively. Glutamate receptors (GluRs) are receptors that bind to glutamate, and they function as ligand-gated ion channels in the central nervous system and mediate transmission in excitatory synapses (41). The subunit composition of AMPARs depends on the conductance properties of Ca2+. Absence of GluR2 subunit promotes permeability to Ca2+, whereas presence of GluR2 inhibits permeability to Ca2+ (42). However, GluR1 and GluR4 subunits also function as Ca2+- permeable AMPARs. Ishiuchi et, al found GluR1 proteins were substantially expressed in most tumor cells, whereas GluR2 was mainly expressed in normal tissues in human glioblastoma samples (43). Furthermore, it has been suggested that blockage of Ca2+ influx through GluR2 expression suppresses migration and induces apoptosis in human glioblastoma cells (44). In addition, knocking down GluR1 inhibits glioma growth (45). Therefore, the conversion of Ca2+- permeable AMPARs to Ca2+- impermeable could be a potential therapeutic target for brain tumors (43). TRIM21 expression is correlated with prognosis, which acts as a tumor suppressor in patients with GBM (46). TRIM21 depletion in GBM enhanced cell proliferation and tumor growth. Lee et, al found that phosphofructokinase 1 (PFK1) expression promotes human glioblastoma progression, while TRIM21 exert anti-tumor effect by mediating poly ubiquitination and degradation of PFK1 (46). Therefore, TRIM21 is a novel target for glioblastoma treatment.

The expression of 24 hub genes in human glioblastomas was validated using single-cell sequencing. Conventional RNA-seq is regularly performed on a bulk level and only measures the average gene expression based on mixed cell populations in samples. Genes that contribute to cell-by-cell variations cannot be detected using conventional RNA-seq data of GBM downloaded from TCGA database (47). However, single-cell RNA-seq (scRNA-seq) profiles for intracellular transcriptome at individual cell level can reveal potential heterogeneous tumors and the composition of glioblastoma tumor microenvironment (48). ScRNA-seq can easily identify highly variable genes in all cell types in the TME of GBM, including the two primary cell types: microglia/macrophages and oligodendrocytes, which are limited in conventional RNA-seq (49). For example, as we have mentioned above, the results of ScRNA-seq revealed that GluR1 and GluR2 were expressed in oligodendrocytes and malignant cells. The expression of Ca2+-permeable GluR confers protection against excitotoxicity and promotes progression of tumor (50). BST2 expression increases in the malignant cells of glioma during tumor progression (51). TLR2 expressed in microglia can promote glioblastoma progression by up-regulating the expression of MT1-MMP in microglia (52). The expression of CCND1 in microglia cells contributes to the differential diagnosis of oligodendrogliomas (53). The use of scRNA-seq to detect the expression of hub genes could significantly help us to accurately understand the function of hub genes in each cell (54). In addition, scRNA-seq demonstrates transcriptional heterogeneity associated with spatial specificity in distinct TME patterns (55). ScRNA-seq has emerged as a revolutionary tool to enhance our understanding of the profiles of hub genes in GBM, and offers insights with implications for both targeted and immune therapies for GBM (49).

In summary, the study identified four types of survival-related immune cells from TCGA database and 24 TME-related hub genes in glioblastoma. The correlation between immune cells and hub genes in patients with GBM was validated using single-cell sequencing data. The results revealed that the hub genes are involved in the development and progression of GBM. Therefore, the candidate genes identified in the study can be used as potential prognostic biomarkers for GBM. However, further studies on the immune cells and hub genes in GBM tumor microenvironment should be conducted to investigate the underlying mechanisms. The present study provides novel insights into the potential association between immune cell TME and GBM prognosis.

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

SH, ZS, JS, and JP designed the study. SH, ZS, XH, TZ, and KH collected and analyzed data. SH, ZS, XH, QZ, and JP wrote the manuscript. All the authors approved the manuscript. All authors contributed to the article and approved the submitted version.


This work was supported by the Natural Science Foundation of Zhejiang Province (project number: LZ20H090002) and the National Natural Science Foundational (project number: 82071285).

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.

Supplementary Material

The Supplementary Material for this article can be found online at:

Supplementary Figure 1 | Overall survival analysis of 22 immune cells based on Kaplan Meier-plotter from the data of CGGA.

Supplementary Figure 2 | Hub genes expression in microglia. Microglia classified from macrophages are marked in dark red, while hub genes are in orange.


1. Touat M, Idbaih A, Sanson M, Ligon KL. Glioblastoma targeted therapy: updated approaches from recent biological insights. Ann Oncol (2017) 28(7):1457–72. doi: 10.1093/annonc/mdx106

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Omuro A, DeAngelis LM. Glioblastoma and other malignant gliomas: a clinical review. Jama (2013) 310(17):1842–50. doi: 10.1001/jama.2013.280319

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Medikonda R, Dunn G, Rahman M, Fecci P, Lim M. A review of glioblastoma immunotherapy. J Neurooncol (2020). doi: 10.1007/s11060-020-03448-1

CrossRef Full Text | Google Scholar

4. Da Ros M, De Gregorio V, Iorio AL, Giunti L, Guidi M, de Martino M, et al. Glioblastoma Chemoresistance: The Double Play by Microenvironment and Blood-Brain Barrier. Int J Mol Sci (2018) 19(10):2879. doi: 10.3390/ijms19102879

CrossRef Full Text | Google Scholar

5. Tomaszewski W, Sanchez-Perez L, Gajewski TF, Sampson JH. Brain Tumor Microenvironment and Host State: Implications for Immunotherapy. Clin Cancer Res (2019) 25(14):4202–10. doi: 10.1158/1078-0432.ccr-18-1627

PubMed Abstract | CrossRef Full Text | Google Scholar

6. 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(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Doncheva NT, Morris JH, Gorodkin J, Jensen LJ. Cytoscape StringApp: Network Analysis and Visualization of Proteomics Data. J Proteome Res (2019) 18(2):623–32. doi: 10.1021/acs.jproteome.8b00702

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Jia D, Li S, Li D, Xue H, Yang D, Liu Y. Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging (Albany NY) (2018) 10(4):592–605. doi: 10.18632/aging.101415

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med (2015) 21(8):938–45. doi: 10.1038/nm.3909

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Cao J, Yang X, Li J, Wu H, Li P, Yao Z, et al. Screening and Identifying Immune-Related Cells and Genes in the Tumor Microenvironment of Bladder Urothelial Carcinoma: Based on TCGA Database and Bioinformatics. Front Oncol (2019) 9:1533:1533. doi: 10.3389/fonc.2019.01533

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Heberle H, Meirelles GV, da Silva FR, Telles GP, Minghim R. InteractiVenn: a web-based tool for the analysis of sets through Venn diagrams. BMC Bioinf (2015) 16(1):169. doi: 10.1186/s12859-015-0611-3

CrossRef Full Text | Google Scholar

12. Huang DW, Sherman BT, Tan Q, Kir J, Liu D, Bryant D, et al. DAVID Bioinformatics Resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res (2007) 35(Web Server issue):W169–75. doi: 10.1093/nar/gkm415

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res (2017) 45(D1):D362–d368. doi: 10.1093/nar/gkw937

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res (2019) 47(D1):D607–d613. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Xu WH, Xu Y, Wang J, Wan FN, Wang HK, Cao DL, et al. Prognostic value and immune infiltration of novel signatures in clear cell renal cell carcinoma microenvironment. Aging (Albany NY) (2019) 11(17):6999–7020. doi: 10.18632/aging.102233

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Chen Q, Yu D, Zhao Y, Qiu J, Xie Y, Tao M. Screening and identification of hub genes in pancreatic cancer by integrated bioinformatics analysis. J Cell Biochem (2019) 120(12):19496–508. doi: 10.1002/jcb.29253

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Wu Y, Zhang S, Yan J. IRF1 association with tumor immune microenvironment and use as a diagnostic biomarker for colorectal cancer recurrence. Oncol Lett (2020) 19(3):1759–70. doi: 10.3892/ol.2020.11289

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res (2017) 77(21):e108–10. doi: 10.1158/0008-5472.can-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Neftel C, Laffy J, Filbin MG, Hara T, Shore ME, Rahme GJ, et al. An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell (2019) 178(4):835–849.e21. doi: 10.1016/j.cell.2019.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Charles NA, Holland EC, Gilbertson R, Glass R, Kettenmann H. The brain tumor microenvironment. Glia (2012) 60(3):502–14. doi: 10.1002/glia.21264

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Chen Z, Hambardzumyan D. Immune Microenvironment in Glioblastoma Subtypes. Front Immunol (2018) 9:1004:1004. doi: 10.3389/fimmu.2018.01004

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Miao X, Leng X, Zhang Q. The Current State of Nanoparticle-Induced Macrophage Polarization and Reprogramming Research. Int J Mol Sci (2017) 18(2):336. doi: 10.3390/ijms18020336

CrossRef Full Text | Google Scholar

23. Komohara Y, Horlad H, Ohnishi K, Fujiwara Y, Bai B, Nakagawa T, et al. Importance of direct macrophage-tumor cell interaction on progression of human glioma. Cancer Sci (2012) 103(12):2165–72. doi: 10.1111/cas.12015

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. Saha D, Martuza RL, Rabkin SD. Macrophage Polarization Contributes to Glioblastoma Eradication by Combination Immunovirotherapy and Immune Checkpoint Blockade. Cancer Cell (2017) 32(2):253–267.e5. doi: 10.1016/j.ccell.2017.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Gordon S. Alternative activation of macrophages. Nat Rev Immunol (2003) 3(1):23–35. doi: 10.1038/nri978

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Levy EM, Roberti MP, Mordoh J. Natural killer cells in human cancer: from biological functions to clinical applications. J BioMed Biotechnol (2011) 2011:676198. doi: 10.1155/2011/676198

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Purdy AK, Campbell KS. Natural killer cells and cancer: regulation by the killer cell Ig-like receptors (KIR). Cancer Biol Ther (2009) 8(23):2211–20. doi: 10.4161/cbt.8.23.10455

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Yang C, Li Y, Yang Y, Chen Z. Overview of Strategies to Improve Therapy against Tumors Using Natural Killer Cell. J Immunol Res (2020) 2020:8459496. doi: 10.1155/2020/8459496

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Golán I, Rodríguez de la Fuente L, Costoya JA. NK Cell-Based Glioblastoma Immunotherapy. Cancers (Basel) (2018) 10(12):522. doi: 10.3390/cancers10120522

CrossRef Full Text | Google Scholar

31. Castriconi R, Daga A, Dondero A, Zona G, Poliani PL, Melotti A, et al. NK cells recognize and kill human glioblastoma cells with stem cell-like properties. J Immunol (2009) 182(6):3530–9. doi: 10.4049/jimmunol.0802845

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Burger MC, Zhang C, Harter PN, Romanski A, Strassheimer F, Senft C, et al. CAR-Engineered NK Cells for the Treatment of Glioblastoma: Turning Innate Effectors Into Precision Tools for Cancer Immunotherapy. Front Immunol (2019) 10:2683:2683. doi: 10.3389/fimmu.2019.02683

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

34. Boix E, Torrent M, Sánchez D, Nogués MV. The antipathogen activities of eosinophil cationic protein. Curr Pharm Biotechnol (2008) 9(3):141–52. doi: 10.2174/138920108784567353

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Curran CS, Evans MD, Bertics PJ. GM-CSF production by glioblastoma cells has a functional role in eosinophil survival, activation, and growth factor production for enhanced tumor cell proliferation. J Immunol (2011) 187(3):1254–63. doi: 10.4049/jimmunol.1001965

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Dustin ML. Integrins and Their Role in Immune Cell Adhesion. Cell (2019) 177(3):499–501. doi: 10.1016/j.cell.2019.03.038

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Harjunpää H, Llort Asens M, Guenther C, Fagerholm SC. Cell Adhesion Molecules and Their Roles and Regulation in the Immune and Tumor Microenvironment. Front Immunol (2019) 10:1078:1078. doi: 10.3389/fimmu.2019.01078

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Serezani CH, Ballinger MN, Aronoff DM, Peters-Golden M. Cyclic AMP: master regulator of innate immune cell function. Am J Respir Cell Mol Biol (2008) 39(2):127–32. doi: 10.1165/rcmb.2008-0091TR

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Takano T, Lin JH, Arcuino G, Gao Q, Yang J, Nedergaard M. Glutamate release promotes growth of malignant gliomas. Nat Med (2001) 7(9):1010–5. doi: 10.1038/nm0901-1010

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Savaskan NE, Seufert S, Hauke J, Tränkle C, Eyüpoglu IY, Hahnen E. Dissection of mitogenic and neurodegenerative actions of cystine and glutamate in malignant gliomas. Oncogene (2011) 30(1):43–53. doi: 10.1038/onc.2010.391

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ozawa S, Kamiya H, Tsuzuki K. Glutamate receptors in the mammalian central nervous system. Prog Neurobiol (1998) 54(5):581–618. doi: 10.1016/s0301-0082(97)00085-3

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Iino M, Goto K, Kakegawa W, Okado H, Sudo M, Ishiuchi S, et al. Glia-synapse interaction through Ca2+-permeable AMPA receptors in Bergmann glia. Science (2001) 292(5518):926–9. doi: 10.1126/science.1058827

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ishiuchi S, Tsuzuki K, Yoshida Y, Yamada N, Hagimura N, Okado H, et al. Blockage of Ca(2+)-permeable AMPA receptors suppresses migration and induces apoptosis in human glioblastoma cells. Nat Med (2002) 8(9):971–8. doi: 10.1038/nm746

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Beretta F, Bassani S, Binda E, Verpelli C, Bello L, Galli R, et al. The GluR2 subunit inhibits proliferation by inactivating Src-MAPK signalling and induces apoptosis by means of caspase 3/6-dependent activation in glioma cells. Eur J Neurosci (2009) 30(1):25–34. doi: 10.1111/j.1460-9568.2009.06804.x

PubMed Abstract | CrossRef Full Text | Google Scholar

45. de Groot JF, Piao Y, Lu L, Fuller GN, Yung WK. Knockdown of GluR1 expression by RNA interference inhibits glioma proliferation. J Neurooncol (2008) 88(2):121–33. doi: 10.1007/s11060-008-9552-2

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Lee JH, Liu R, Li J, Zhang C, Wang Y, Cai Q, et al. Stabilization of phosphofructokinase 1 platelet isoform by AKT promotes tumorigenesis. Nat Commun (2017) 8(1):949. doi: 10.1038/s41467-017-00906-9

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet (2009) 10(1):57–63. doi: 10.1038/nrg2484

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science (2014) 344(6190):1396–401. doi: 10.1126/science.1254257

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Tirosh I, Suvà ML. Dissecting human gliomas by single-cell RNA sequencing. Neuro Oncol (2018) 20(1):37–43. doi: 10.1093/neuonc/nox126

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Yoshioka A, Ikegaki N, Williams M, Pleasure D. Expression of N-methyl-D-aspartate (NMDA) and non-NMDA glutamate receptor genes in neuroblastoma, medulloblastoma, and other cells lines. J Neurosci Res (1996) 46(2):164–78. doi: 10.1002/(sici)1097-4547(19961015)46:2<164::aid-jnr4>;2-f

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Wainwright DA, Balyasnikova IV, Han Y, Lesniak MS. The expression of BST2 in human and experimental mouse brain tumors. Exp Mol Pathol (2011) 91(1):440–6. doi: 10.1016/j.yexmp.2011.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Vinnakota K, Hu F, Ku MC, Georgieva PB, Szulzewsky F, Pohlmann A, et al. Toll-like receptor 2 mediates microglia/brain macrophage MT1-MMP expression and glioma expansion. Neuro Oncol (2013) 15(11):1457–68. doi: 10.1093/neuonc/not115

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Bosone I, Cavalla P, Chiadò-Piat L, Vito ND, Schiffer D. Cyclin D1 expression in normal oligodendroglia and microglia cells: its use in the differential diagnosis of oligodendrogliomas. Neuropathology (2001) 21(3):155–61. doi: 10.1046/j.1440-1789.2001.00389.x

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Yip SH, Sham PC, Wang J. Evaluation of tools for highly variable gene discovery from single-cell RNA-seq data. Brief Bioinform (2019) 20(4):1583–9. doi: 10.1093/bib/bby011

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Venteicher AS, Tirosh I, Hebert C, Yizhak K, Neftel C, Filbin MG, et al. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science (2017) 355(6332):eaai8478. doi: 10.1126/science.aai8478

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: glioblastoma, tumor microenvironment, immune infiltration, immune therapy, TCGA

Citation: Huang S, Song Z, Zhang T, He X, Huang K, Zhang Q, Shen J and Pan J (2020) Identification of Immune Cell Infiltration and Immune-Related Genes in the Tumor Microenvironment of Glioblastomas. Front. Immunol. 11:585034. doi: 10.3389/fimmu.2020.585034

Received: 19 July 2020; Accepted: 28 September 2020;
Published: 20 October 2020.

Edited by:

Xiaoxing Xiong, Renmin Hospital of Wuhan University, China

Reviewed by:

Rongjing Ge, Bengbu Medical College, China
Xianli Lv, Tsinghua University, China

Copyright © 2020 Huang, Song, Zhang, He, Huang, Zhang, Shen and Pan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Jian Shen,; Jianwei Pan,

These authors have contributed equally to this work