Investigation of Genetic Determinants of Glioma Immune Phenotype by Integrative Immunogenomic Scale Analysis

The immunosuppressive mechanisms of the surrounding microenvironment and distinct immunogenomic features in glioblastoma (GBM) have not been elucidated to date. To fill this gap, useful data were extracted from The Cancer Genome Atlas (TCGA), the Chinese Glioma Genome Atlas (CGGA), GSE16011, GSE43378, GSE23806, and GSE12907. With the ssGSEA method and the ESTIMATE and CIBERSORT algorithms, four microenvironmental signatures were used to identify glioma microenvironment genes, and the samples were reasonably classified into three immune phenotypes. The molecular and clinical features of these phenotypes were characterized via key gene set expression, tumor mutation burden, fraction of immune cell infiltration, and functional enrichment. Exhausted CD8+ T cell (GET) signature construction with the predictive response to commonly used antitumor drugs and peritumoral edema assisted in further characterizing the immune phenotype features. A total of 2,466 glioma samples with gene expression profiles were enrolled. Tumor purity, ESTIMATE, and immune and stromal scores served as the 4 microenvironment signatures used to classify gliomas into immune-high, immune-middle and immune-low groups, which had distinct immune heterogeneity and clinicopathological characteristics. The immune-H phenotype had higher expression of four immune signatures; however, most checkpoint molecules exhibited poor survival. Enriched pathways among the subtypes were related to immunity. The GET score was similar among the three phenotypes, while immune-L was more sensitive to bortezomib, cisplatin, docetaxel, lapatinib, and rapamycin prescriptions and displayed mild peritumor edema. The three novel immune phenotypes with distinct immunogenetic features could have utility for understanding glioma microenvironment regulation and determining prognosis. These results contribute to classifying glioma subtypes, remodeling the immunosuppressive microenvironment and informing novel cancer immunotherapy in the era of precision immuno-oncology.


INTRODUCTION
Gliomas are the most common and malignant primary tumors in the central nervous system (CNS) and have a highly invasive nature (1,2). The discovery of the lymphatic system in the CNS has aroused inspiration to provide a novel theoretical foundation and new prospects for immunotherapy in brain tumors, and previous literature has demonstrated the mutual interactions between gliomas and the immune system (3,4). Multiple related biological processes influencing CNS immune surveillance, such as the PI3K/ Akt pathway, FAK, the IGF pathway, the STAT3 pathway, chemokines, HIF-1a, IL-6, TGF-b, PD-1/PD-L1, and CTLA-4, could individually impact immunosurveillance (5)(6)(7)(8). Since entering the era of precision oncology, the molecular profiles of gliomas have been well studied. Mutations in the isocitrate dehydrogenase 1 (IDH1) gene, 1p/19q codeletion, methylguanine methyltransferase (MGMT) promoter methylation, tumor protein 53 (TP53), and telomerase reverse transcriptase (TERT) promoters are becoming treatment targets or prognostic biomarkers (9)(10)(11). Monoclonal antibodies (mAbs) against PD-1/PD-L1 show satisfying overall survival (OS) in melanoma and non-small cell lung cancer (NSCLC), but there is limited survival benefit in glioma (12). The unique immune-privileged microenvironment due to the inherent expression of immunosuppressive cytokines, such as PD-1, TGF-b, and IL-10, and the lack of antigen-presenting cells (APCs) in the CNS present obstacles for the efficacy of immunotherapy in glioblastoma (GBM) (13). The development of more novel and effective therapies will require a deep understanding of the tumor immunosuppressive microenvironment.
Direct interactions between tumor and immune cells can result in suppression of natural killer (NK) cell activity mediated by HLA molecules (including HLA-E and HLA-G) (14), immune cell apoptosis via tumor necrosis factor receptor superfamily member 6 (TNFRSF6, known as FAS) and the FAS-ligand interaction (15), or triggering of inhibitory T cell checkpoints by PD-L1 (16). The hypofunctional state of T cells known as T-cell exhaustion was identified by the accumulation of coinhibitory checkpoints (17). Of note, the paucity of T cells in the glioma microenvironment is striking in contrast to the levels in other "hot tumors", and some studies have suggested that glioma-associated myeloid cells are immunosuppressive with an M2-like phenotype and require peripheral dendritic cells (DCs) to elicit an immune response (18). Indeed, the exact mechanism of immune suppression is still obscure. In this study, we employed 2,466 samples to properly classify glioma into immune phenotypes according to distinct immunogenomic features based on microenvironment-related genes. Then, we validated and identified microenvironmental processes, explored immune alterations, and characterized immunosuppressive mechanisms. The immune landscape may inspire glioma subtype classification, remodeling of the immunosuppressive microenvironment and development of new therapies.

Immune Phenotype Classification
In the glioma microenvironment, immune and stromal cells are two key types of nontumor components and have been indicated to be significant for the diagnosis and prognosis of tumors. Yoshihara et al (19) designed the ESTIMATE algorithm to compute immune and stromal cell scores to predict the infiltration of these nontumor cells. The authors used ESTIMATE to evaluate immune scores, ESTIMATE scores, stromal scores and tumor purity scores in each tumor sample with the aim of determining the immune infiltration level.
Single-sample gene set enrichment analysis (ssGSEA), which assisted in quantifying the enrichment level of an immune cell/ signature, pathway or biological process in a tumor sample, was used to assess the gene score of every gene set for every sample (20). The enrichment-related score represented the level at which the genes in the gene set were synchronously up-or downregulated in the sample. The infiltration of immune cells in the microenvironment was determined by 29 immune cell types: NK cells, effector memory CD4+ T cells, activated B cells, monocytes, memory B cells, activated CD4+ T cells, type 2 T helper cells, dendritic cells, neutrophils, macrophages, effector memory CD8+ T cells, myeloid-derived suppressor cells (MDSCs), immature B cells, mast cells, and regulatory T cells, and glioma samples were hierarchically clustered into "immune-high (immune-H)", "immune-middle (immune-M)" and "immune-low (immune-L)" groups. Separation of gene expression patterns between immune phenotypes was evaluated by the principal component analysis (PCA) algorithm with the PCA1, PCA2, and PCA3 top three dimensions (21). Visualization was performed with the "GSVA", "GSEABase", "ComplexHeatmap", "estimate", and "ggplot" public packages.

Quantification of Molecular and Genomic Features
Tumor mutation burden (TMB) was defined as the total count of somatic mutations per megabase in each tumor sample. We used the MATH algorithm (22), which assessed the width of the allele frequency distribution, to evaluate the intratumor heterogeneity level of tumor samples. Further intratumor heterogeneity scores were quantified using the function "math. Score" in the "maftools" package with downloaded "maf" files based on the hg19 sequencing platform. Comparisons of the somatic mutations and SNP sites among immune phenotypes in distinct populations (low-grade glioma (LGG) and GBM samples) were displayed to investigate the discrepancies by the "maftools" and "GenVisR" packages.

Survival Analysis
Available clinicopathological factors (e.g., sex, age, treatment options, histological subtype, and classic mutations) were collected from the TCGA and CGGA datasets to estimate the association between these factors as well as the immune phenotypes and prognosis with univariable and multivariable Cox analysis (uniCox, multiCox) and proportional hazard models. We compared survival differences among immunespecific phenotypes of glioma in distinct groups using Kaplan-Meier curves and the log-rank test with normalized clinical data.

Estimation of the Proportions of Immune Cell Types
CIBERSORT is an algorithm designed to characterize the cell composition of complex tissues based on their gene expression profiles, and it is highly consistent with real-life estimations in many cancers. A leukocyte gene signature matrix employing 547 genes, which was defined as LM22, was used to quantify 22 immune cell types (23). These 22 types of immune cells mainly include myeloid subtypes, NK cells, plasma cells, naive and memory B cells and T cells. We used the CIBERSORT method to investigate the fraction of the 22 immune cell types in each derived phenotype and identify the characteristics of infiltrating cells in the glioma microenvironment.
Identification of a Gene Signature for Exhausted CD8+ T Cells CD8+ T lymphocytes are regarded as a critical component of antitumor immunity, while immune invasion often occurs during the development of T cell exhaustion, characterized by the progressive accumulation of coinhibitory checkpoints, including PD-1, PD-L1, CTLA-4, TIM-3, and LAG-3 (17). We defined a gene expression signature of exhausted CD8+ T cells with integrative bioinformatics through publicly available NSCLC data considering the data quality and availability. We obtained an RNAseq dataset of intratumoral CD8+ T cells showing high or no PD-1 (PDCD1) expression in a published study (24), and we generated an upregulated PD-1-positive gene list from another previous study (25). Pearson correlation analysis was conducted using the upregulated PD-1-positive gene list in the TCGA (microarray+ RNA-seq cohort) and CGGA (microarray+ RNA-seq cohort) datasets with an adjusted P-value < 0.05 and |correlation efficiency| > 0.25 as the eligibility criteria. In total, a 5-gene signature was identified in the glioma database, and an exhausted CD8+ T cell (GET) score was quantified in a tumor by conducting ssGSEA to obtain the ssGSEA score. In combination with clinical and molecular profiles, the prognostic and predictive values of the GET score were determined through different immune phenotypes.

Correlation and Functional Analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed on genes differentially expressed between the immune-high and immunelow groups. Gene set enrichment analysis (GSEA) was carried out to identify the group of genes enriched either in the immune-high or immune-low group with cutoffs of a P-value < 0.1 and a false discovery rate (FDR) < 0.05 (26). Gene set variation analysis (GSVA) is a nonparametric and unsupervised method estimating the variations of samples in analyzed datasets in pathways and biological process (27). The gene sets of "c2.cp.kegg.v6.2.-symbols" used were captured from the MSigDB website for GSVA, with an adjusted P-value < 0.05 considered statistically significant. Correlation plots were constructed that primarily focused on the interactions between IDH1 and other key immune-related genes identified from the GSEA with a P filter = 0.001. A Sankey diagram was constructed to show the correlations between checkpoints and the GSEA-identified genes. Visualization of the unction analyses was realized via the "circlize" (28), "circus" (29), "clusterProfiler", and "ggalluvial" (30) packages.

Prediction of the Chemo/Targeted Therapy Response
Intended chemotherapeutic and targeted responses of glioma samples were evaluated by the largest publicly available pharmacogenomics database (Pharmaceutical Sensitivity Genomics in Cancer (GDSC) https://www.cancerrxgene.org/) (31). GDSC contains drug sensitivity information from nearly 75000 experiments and responses to 138 anticancer drugs across almost 700 cell lines. The database provides a unique source relevant to mainstream drug sensitivity and genomic datasets to inspire new discoveries on cancer therapeutic biomarkers. GDSC is also utilized due to its visualization capability. The evaluation procedure was conducted via the R software package "pRRophetic", half-maximal inhibitor concentration (IC50), and the evaluation accuracy was determined by ridge regression and 10-fold cross-validation using the GDSC dataset (32,33). Different chemotherapeutic and targeted responses among the three phenotypes were analyzed by one-way analysis of variance (ANOVA) or the Kruskal-Wallis test (K-W test) based on the results of the normal distribution criteria test. The response to commonly used chemotherapy or targeted therapies was compared according to immune phenotype, although some drugs were not formally approved for utility in brain tumors.

Peri-Tumoral Edema Characteristics
To detect the variations in some radiomics features of classified immune phenotypes, MR images (MRIs) of patients from the TCGA dataset were obtained from the Cancer Imaging Archive. TCGA-GBM and TCGA-LGG cohorts in the Cancer Imaging Archive (http://www.cancerimagingarchive.net) were specifically selected and matched with existing results. Eligible tumor contrast enhancement images were determined after a discussion with three skilled neurosurgeons (Zhao B, Xing H, Wang Y) on the author list. Radiomics features of tumors included tumor size, enhancement, noncontrast-enhancing tumor (nCET), necrosis, edema, cysts, multifocality, contact with ventricles or neocortex and location based on a previous study (34). Features such as multifocality, enhancement, location and edema were revealed to have molecular signature correlations with glioma, such as IDH mutation or MGMT promoter methylation; edema and necrosis were regarded as poor survival imaging markers (34,35). Edema associated with both molecular phenotypes and prognosis was the focus of investigation to facilitate identification of noninvasive acquired markers and features of the classified glioma phenotypes. A mild (or no) region of edema (-) was regarded as edema extending ≤ 1 cm from the margin of the tumor; otherwise, it was treated as moderate to severe (+) (36). The evaluations were all based on eligible T2-weighted images.

Statistical and Bioinformatics Analyses
Statistical analyses were conducted using R software (version 3.5.3), and other statistical methods are mentioned throughout the article. Bioinformatics analysis was conducted mainly following the methods of Thorsson et al (37). A two-sided P < 0.05 was considered to be significant unless otherwise specified. The public packages used are mentioned throughout this paper.

ssGSEA and Independent Immune Phenotype Classification
After excluding the normal tissues (5 normal samples in the TCGA RNA-seq database), tumor samples with distinct extension of inflammatory cell infiltration were classified into "immune-L", "immune-M" and "immune-H" phenotypes with ssGSEA incorporating 29 types of immune cell lineages, such as helper T cells, cytotoxic T cells, myeloid cells, monocytes, NK cells, dendritic cells, and T cells.

Each Phenotype Has Distinct Immunogenetic Features
Four immune scores were employed. From the ESTIMATE algorithm, the immune-H phenotype was revealed to have a higher ESTIMATE score, immune score and stromal score and a lower tumor purity score than the immune-M and immune-L phenotypes. Statistical comparisons showed that there were significant differences between the immune-H and immune-L phenotypes (Wilcoxon P-value < 0.001) related to these immune signatures ( Figure 2).
Checkpoint biomarkers are involved in tumor subtype classification, prognosis prediction and immunotherapy therapy  response evaluation. We observed that most checkpoints were differentially expressed. Such biomarkers were more highly expressed in the immune-H phenotype than in the immune-M and immune-L groups. CD200 was highly expressed in the immune-L phenotype (K-W test P value < 0.001) ( Figure 3A). HLA genes took important roles in innate immunity and tumor immune microenvironment regulation, these family genes had significantly different expression among phenotypes, with the immune-H group exhibiting significantly higher expression than the other two groups ( Figure 3B). Furthermore, the immune-L showed higher expression of TP53, EGFR, NF1, PDGFRA, and RB1than immune-H phenotype, which suggested the converged axes of P53, tumor suppressive Rb and MAPK/PI3K were potentially activated in immune-L phenotype. IDH-mutant glioma with ATRX and TERT mutations was always associated with favorable survival ( Figure 3C). Good separation between the immune-H and immune-L phenotypes was confirmed by PCA ( Figures 3D-F). Based on the above results, the immune-H phenotype may be more sensitive to classic checkpoint immunotherapy than the others, while the immune-L phenotype was associated with longer survival and better prognosis. P-value = 2.056e-12), the immune-H phenotype exhibited unfavorable survival compared with the immune-L phenotype ( Figures 4E, F). Similar findings were consistent and confirmed in the TCGA RNA-seq (P-value < 0.001), CGGA microarray (P-value = 1.135e-5) and CGGA RNA-seq cohorts (P-value = 8.882e-16), but the results were not significant in the TCGA microarray cohort due to limitations derived from the sample number (P-value = 0.186) ( Figures 4G-J). For subgroup analyses conducted in the CGGA cohort, in the LGG and primary glioma patients, there were significant survival differences between the immune-H, -M and -L phenotypes (log-rank P-value = 7.346e-4; P-value = 9.448e-14, respectively). The prognostic value was not obvious in the GBM (P-value = 0.928) or recurrent subpopulations (P-value = 0.658) ( Figures 4K-N). These results were contrary to those of previous studies on other cancer types, including breast cancer (38), gastric cancer (39) and head and neck squamous cell cancer (40), which indicated the specificity of the association between tumor immunity and clinical outcomes in glioma, the microenvironment of which is regarded as rather immunosuppressive and refractory. Additionally, intrinsic limitations associated with sample size and variation of ethnicity among the used databases or cohorts should be acknowledged.

Infiltrating Immune Cell Fractions and Correlations
Through the CIBERSORT algorithm, M2, M1, and M0 macrophages, monocytes, DCs, and subsets of B and T cells (CD4+ and CD8+) were distinguished in the glioma microenvironment ( Figures 5A, B). The results derived from ESTIMATE and CIBERSORT classified the glioma samples into three immune phenotypes, which had similar characteristics to those of the previously identified phenotypes. Correlations between each type of immune cell illustrated that the most negative correlations were found among M0 macrophages, monocytes, M2 macrophages, DCs (activated and resting) and helper T cells. These results suggested that myeloid cells highly  participated in the immunosuppressive glioma microenvironment ( Figures 5C, D). Comparing the proportion of infiltrating immune cells, the immune-H phenotype was revealed to have higher proportions of all analyzed immune cells, except M2 macrophages, activated mast cells, monocytes, neutrophils and resting memory CD4+ T cells ( Figure 5E).

Construction of the Exhausted CD8+ T Cell Signature
Exhausted CD8+ T cell levels were recognized to be uniquely regulated by distinct PD-1 upregulation. With transcriptional profiles of CD8+ T lymphocytes and upregulated PD-1-positive genes captured from previous studies (24,25), correlation analyses were carried out in the involved datasets, in which five genes meeting the selection criteria were selected and termed GET signature. The GET signature included PDCD1, CD27, ICOS, RUNX2, and CXCR6, which are closely linked to T cell dysfunction and coregulation ( Figure 6F). The GET score of each tumor sample was established with the ssGSEA method. To quantitatively illustrate the status of exhausted CD8+ T cells in each immune phenotype, we compared the distribution of the GET score in different phenotypes. We did not observe significant differences in the GET score between the immune-L, -M and -H phenotypes ( Figure 6A, Supplementary Online File 5). Correlations between the defined GET score and immune score, ESTIMATE score, stromal score and tumor purity were assessed, and no tight correlation was found among these signatures ( Figures 6B-E). The results from the TCGA microarray dataset seemed to vary slightly from the results in other datasets, and the lack of CD8+ T cells in the glioma microenvironment and the failure of immune surveillance against tumor cells were likely causes of these effects (41). Patients with a higher GET score in the total CGGA cohort had a more favorable prognosis than those with a lower GET score (HR: 1.38, 1.20-1.60; P-value = 1.25e-5), and the results were not significant in the total TCGA cohort ( Figures 6G, H) (Supplementary Online File 6). Confirmatively, nearly all of the constructed GET signatures were mainly related to inflammatory components, lymphocyte functions and immune cell signaling ( Figure 6I). To date, crosstalk between the GET signature and other molecular and clinicopathological factors is being warmly discussed in glioma, and more evidence is needed in the future.

Functional Enrichment Analysis of Phenotype-Associated Genes
In subsequent functional analyses of the biological processes of the identified microenvironment-related genes in the immune phenotypes, metagenes chosen as classifier gene sets for the immune-H over the immune-L phenotype in GSEA were significantly enriched in immune-related GO terms such as dendritic cell antigen processing and presentation, immunoglobulin processes, regulation of T cell chemotaxis, and T helper cell lineage (P-value and Benjamini P-value < 0.05); metagenes were significantly enriched in pathways related to immune-related graft-versus-host disease, the hematopoietic cell lineage, and the IL-17 signaling pathway (P-value and Benjamini P-value < 0.05) according to pathway GSEA ( Figures 7A-D). Bubble plots can be found in Appendix Figure A1. The cluster maps display whole gene  (Figures 7E-H). GSVA showed enrichment discrepancies in several immune-related pathways, including antigen processing, primary immunodeficiency, the B/T cell receptor signaling pathway, NK cell cytotoxicity, and leukocyte transendothelial migration (Figures 7I, J). The Sankey diagram shows the links between checkpoint molecules and their correlated genes in glioma ( Figure 7K).

Genomic Alterations, Tumor Mutation Burden, and Histological Characteristics
Compared with other immune phenotypes, immune-L was found to have a higher proportion of IDH-mutant patients ( Figures 8A, D); the immune-H phenotype seemed to have a higher proportion of recurrent glioma but a lower rate of primary patients ( Figures 8B, E); more GBM samples existed in immune-H, and more LGG was associated with the immune-L phenotype ( Figures 8C, F). Detailed data are presented in Table 2. We analyzed the distribution differences of somatic mutations and SNPs among the immune phenotypes using data from the TCGA project. Figures 9A, B displays recurrent SNP sites (N > 5) in chromosome models in LGG and GBM. Sites marked by orange and red are high-mutant SNP sites, while those marked by navy and green are low-mutant SNP sites. Major mutant genes and mutation types were different among immune phenotypes in combination with glioma grade (Figures 9C, D). In addition, GBM presented more extensive TMB than LGG, with details in the left bar plots and scatter plots in Figure 9E.

Phenotypes Predicting Response to Antitumor Drugs and Peritumor Edema
Chemotherapy and targeted therapy are both standard treatments for glioma. The response to commonly used antitumor drugs was evaluated among three immune phenotypes. The expected model using the GDSC dataset was trained by ridge regression, and the level of prediction accuracy was evaluated by 10-fold crossvalidation. The treatment-related IC50 for each tumor sample in TCGA was properly estimated based on a predictive model of these drugs. There were significant differences in the response to several drugs, and the immune-L phenotype was more sensitive to bortezomib (K-W P < 2.2e-16), cisplatin (P = 5.3e-15), docetaxel (P < 2.2e-16), lapatinib (P < 2.2e-16), and rapamycin (P = 3.3e-8); however, the immune-H phenotype was more sensitive to paclitaxel (P = 3.1e-10) and sorafenib (P = 0.0053) ( Figure 10A).
As a marker of inflammation, edema is a common pathophysiological entity surrounding gliomas. Herein, we compared the edema differences between the immune-L and immune-H phenotypes to assess the correlations. It was noted that immune-H phenotype gliomas displayed more severe edema than immune-L phenotype gliomas ( Figure 10B). The present results suggested that peritumoral edema is also a probable marker to reflect the variations between immune phenotypes. The analysis process used in this study is shown as a flow chart in Figure 11.

DISCUSSION
Immunotherapy has been confirmed to be effective in other types of cancers except glioma, as glioma features a relatively immuneprivileged microenvironment. With the aim of elucidating the immunosuppressive mechanism, in this research, we enrolled 2,466 glioma samples from 6 datasets and classified these samples into 3 immune phenotypes with distinct immunogenetic features. The immune-H phenotype has higher immune cell lineage infiltration and higher ESTIMATE, immune and stromal scores than the immune-L and immune-M phenotypes. Most HLA family genes and checkpoint molecules were significantly highly expressed in the immune-H phenotype; otherwise, some specific genes were highly expressed. Overall, patients with the immune-H phenotype will have a poor prognosis compared with those with the immune-L phenotypes, but this result was limited by the sample size. A five-gene GET signature including PDCD1, CD27, ICOS, RUNX2, and CXCR6 was established, and no significant differences in the GET score between immune phenotypes were observed. Patients with a high GET score seemed to have a better prognosis. A response difference was noticed among the phenotypes to several antitumor drugs. Immune-H was observed to have more severe peritumor edema than immune-L in representative T2 images. Survival differences among the classified immune phenotypes of glioma were in contrast to those of some other cancer types reported previously, such as triple-negative breast cancer (38), gastric cancer (39), and head and neck squamous cancer (40). A potential reason is that the inflammatory microenvironment upregulated the tumor progressive nature and deteriorated glioma invasion and development (41,42). The success of immunomodulatory therapy is widespread among diverse cancer types, which stimulates our interest in characterizing TME immune cell infiltration in glioma. The immune-H phenotype may be involved in immunosuppressive activities, including immunosuppressive checkpoints ( Table 3), expression of tumor-supportive macrophage chemotactic and polarizing molecules and immune-suppressive pathway signaling (the IL-10 signaling pathway). The IL-10 pathway downregulates DC activation and IL-12 production and inhibits the cytotoxic T cell response during chemotherapy. Macrophage activation is also suppressed by IL-10 to inhibit the immune response (43). Importantly, there is large heterogeneity in the TME of different glioma genetic subtypes, and enriched tumor-associated macrophages (TAMs) participate in the promotion of glioma invasion, angiogenesis, tumor metastasis and immune suppression through intracellular and extracellular mediators (44). Glioma with IDH mutation status was shown to have low levels of infiltrating T cells and a higher ratio of TAMs derived from microglia (45). Although TAMs have distinct genetic profiles involving canonical M1 (antitumorigenic) and canonical M2 (protumorigenic) polarization, they show increased anti-PD-1 resistance-associated genes and predict poor survival (46,47). Additionally, immunosuppressive chemokines/cytokines in the TME released by the tumor itself, such as through the TGF-b pathway, also block antitumor immunity activation (48). TIM-3 (T cell immunoglobulin mucin receptor 3) has an immunosuppressive effect in glioma, which may be due to the unique presence of TIM-3+ Tregs in tumor tissue (49). Furthermore, TIM-3 does not contain any immunoreceptor tyrosine-based inhibition motifs (ITIMs), which are necessary for avoiding major deficiencies in immunotherapy (50). VISTA (V-type immunoglobulin domain-containing suppressor of T cell activation) is a newly found checkpoint that restricts T cell activation by shaping the naive CD4+ T cell compartment (51). Therapeutics targeting VISTA curb the development of graftversus-host disease and promote the death of naive CD4+ T cells; thus, VISTA can be regarded as a distinctive immunotherapy molecule (51,52). Indeed, growing evidence suggests that dysfunctional CD8+ T cells incorporate heterogeneous subpopulations such as progenitor and terminally exhausted cells, and discrete functions in immunotherapy or the microenvironment need to be better elucidated (53). Clinical trials regarding Checkmate 143 (NCT02017717), Checkmate 498 (NCT02617589), and Checkmate 548 (NCT02667587) did not suggest a profound survival benefit from immunotherapy in glioma/GBM, with only some clinical advantages reported in some case reports; indeed, GBM typically has a relatively low mutational load and a paucity of T cell infiltration compared with other cancers (12,54).
Similar to other studies, Chen and his colleagues (55) used ssGSEA to identify the immune microenvironment of glioma, and they did not classify glioma samples into immune phenotypes or detect the corresponding microenvironmental features of the phenotypes; however, they detected eight glioma microenvironment-associated genes, CCDC109B, EMP3, ANXA2, Chi-square test was conducted to compare these differences between immune phenotypes. IDH MT, IDH Mutant; IDH WT, IDH Wild Type; LGG, low grade glioma; GBM, glioblastoma. CLIC1, TIMP1, VIM, LGALS1, and RBMS1, and constructed a prognostic model with them through integrative omics data points. They validated the immunosuppression of LGALS1 in in vitro experiments. Our findings based on large genomic data help characterize the glioma microenvironment and understand tumor immune complexity. The ESTIMATE, immune, stromal, and tumor purity scores can be used properly in both basic and translational medicine to help identify glioma subtypes. Work investigating the immunosuppressive mechanisms of glioma implies that microenvironments lacking T cells feature immunosuppressive biological processes carried out by a series of immune cells; more knowledge of immune cell infiltration will inform strategies to remodel the immunosuppressive microenvironment and will aid the identification of more therapeutic targets.
Patients with the immune-H phenotype were more prone to developing a poor prognosis compared with others; thus, we may properly predict the prognosis of glioma patients with immune phenotypes. Our findings also suggest that immunotherapy will be effective in immune-H patients, who are more sensitive to checkpoint-related immunotherapy (56). Recent evidence showed that samples with high TMB could exhibit a durable response to PD-1/PD-L1 immunotherapy (57), and current findings indirectly confirmed the value of TMB in predicting immunotherapeutic outcomes of established immune phenotypes. Translational research indicated that a high TMB status may yield a long-term response and durable survival benefit (58). The presented results provide a novel perspective on immune signatures in the genetic TMB, the microenvironment and roles in immune checkpoint The immune-L phenotype was more sensitive to bortezomib (P < 2.2e-16), cisplatin (P = 5.3e-15), docetaxel (P < 2.2e-16), lapatinib (P < 2.2e-16), rapamycin (P = 3.3e-8); the immune-H phenotype was more sensitive to paclitaxel (P = 3.1e-10) and sorafenib (P = 0.0053). (B) Representative images of the differences in the extent of peri-tumoral edema in TCGA cohort patients. Immune-H phenotype significantly possessed more-severe edema than immune-L.  blockade treatment and inspired the exploration of fresh neoepitopes. Immune phenotype classification highlights the importance of individualized treatments and provides potential methods to be used in further clinical trials related to glioma immunotherapy. We believe that with the current Pan-Cancer Analysis of Whole Genomes (PCAWG) project involving classic glioma microenvironment biomarkers (i.e., IDH1), researchers will identify more specialized features of cancer immune genomes (59).

CONCLUSIONS
Glioma samples can potentially be classified into "immune-H", "immune-M" and "immune-L" phenotypes, which exhibit distinct immunogenetic features. The immune-H phenotype is associated with higher ESTIMATE, immune and stromal scores but poorer survival than the immune-L phenotype. HLA and checkpoint family genes are relatively highly expressed in patients with the A summary of the ligands, immune-related expression pattern, biological function, and molecular mechanisms is reviewed for selected costimulatory and coinhibitory receptors. Molecular functions (i.e., downstream signaling) reflect predominant currently known mechanisms, but additional mechanisms are likely to contribute significantly. NK, natural killer; NKT, natural killer T cell; TFH, T follicular helper; TRAF, tumor necrosis factor receptor-associated factors; DC, dendritic cell.
immune-H phenotype. The GET signature cannot effectively reveal the discrepancies among immune phenotypes, and aggressive peritumor edema was displayed in immune-H compared with immune-L phenotypes. Our immunogenetic pipeline characterizes the glioma microenvironment and properly identifies patients who are more sensitive to chemo/targeted therapy and are likely to have better survival. These results possibly facilitate new therapeutic development and advance precision oncology, limited by the observational nature, the experimental profile should be highlighted in the future.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
All authors designed and conducted this review. BZ, YKW, and YaW, wrote the paper. YKW helped the study design. CD and YuW revised the statistical methodology. YuW and WM had primary responsibility for the final content. All authors contributed to the article and approved the submitted version. Notably, YuW and WM equally share the corresponding authorship.