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

Abstract

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-1α, IL-6, TGF-β, PD-1/PD-L1, and CTLA-4, could individually impact immunosurveillance (5–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–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-β, 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.

Methods

Data Acquisition and Filtration

Data from glioma patients from six mRNA databases were extracted from TCGA database (RNA-sequencing (RNA-seq) for GBM, n = 169, microarray, n = 539) (http://cancergenome.nih.gov/), the CGGA database (RNA-seq, n = 1018, microarray, n = 301) (http://www.cgga.org.cn), the GSE16011 database (microarray, n = 276), the GSE43378 database (microarray, n = 50), the GSE23806 database (microarray, n = 92) and the GSE12907 database (microarray, n = 21). Complete clinical information was obtained from TCGA (http://cancergenome.nih.gov/, n = 708) and GCCA (http://www.cgga.org.cn, n = 1319). Somatic mutations and single nucleotide polymorphisms (SNPs) of gliomas were obtained from the TCGA database (http://cancergenome.nih.gov/, n=901, gene number n = 13,389). RNA-seq data downloaded in FPKM values from TCGA were normalized and transformed into transcripts per kilobase million values. RNA expression of gliomas was assessed with the Affymetrix microarray platform in the Gene Expression Omnibus (GEO) database (GSE16011, GSE43378, GSE23806, and GSE12907). Data were filtered to exclude samples without mRNA expression or clear histology, and the genomic data were normalized and analyzed within lanes, between lanes, and per quantile using the “limma” and “DESeq2” R packages. In this study, TCGA and CGGA were mainly treated as the training sets, and GEO databases were regarded as the validation sets.

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 immune-specific 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 RNA-seq 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 immune-low 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/TargetedTherapy 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.

Results

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. The numbers of samples falling into the immune-L, immune-M, and immune-H phenotypes were 283, 234 and 21 in the TCGA microarray data; 129, 8 and 32 in the TCGA-GBM RNA-seq data; 105, 90 and 106 in the CGGA microarray cohort; 413, 425 and 180 in the CGGA RNA-seq cohort; 112, 162 and 2 in GSE16011; 28, 16 and 6 in GSE43378; 87, 2 and 3 in GSE23806; and 9, 10 and 2 in GSE12907, respectively (Figure 1).

Figure 1

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).

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.

Figure 3

The Immune-H Phenotype Is Associated With a Poor Prognosis

Clinical and molecular features of the immune-specific phenotypes of glioma are displayed in complex heatmaps (Figures 4A–D, Supplementary Online Files 1–4). Treatment options and histological characteristics seemed to have more prognostic influence, and patients who had received corresponding chemotherapy (including adjuvant temozolomide (TMZ) therapy) or radiotherapy or who had a lower tumor grade and malignancy were observed to have favorable survival. The results are summarized in standardized Table 1. In the TCGA (n = 701) (log-rank P-value = 0.031) and CGGA cohorts (n = 1281) (log-rank 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.

Figure 4

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).

Figure 5

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.

Figure 6

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 clusters and enriched GO terms, and the GO chord plots show relevant GO terms for the classic PD1/PDCD1, CTLA-4, TIGIT, VISTA/VISR, and LAG-3 molecules (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).

Figure 7

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. Surprisingly, no obviously significant discrepancies in TMB were found between the immune-H and immune-L phenotypes in the TCGA microarray cohort (P = 0.047, median log₂(TMB), 0.385 vs. 0.464) and TCGA RNA-seq cohort (P = 0.100, median log₂(TMB), 0.357 vs. 0.447) (Figures 8G, H, Supplementary Online File 7).

Figure 8

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.

Figure 9

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 cross-validation. 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).

Figure 10

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.

Figure 11

Discussion

Immunotherapy has been confirmed to be effective in other types of cancers except glioma, as glioma features a relatively immune-privileged 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-β 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 graft-versus-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, 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 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 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.

Funding

The work is supported by Chinese Academy of Medical Sciences Innovation Fund for Medical Sciences, number of grants (2016-I2M-2-001) and Beijing Municipal Natural Science Foundation, number of grants (7202150, 19JCZDJC64200(Z)).

References

• 1

  SiegelRLMillerKDJemalA. Cancer Statistics, 2020. CA Cancer J Clin (2020) 70(1):7–30. doi: 10.3322/caac.21590

  • CrossRef
  • Google Scholar

• 2

  JiangTMaoYMaWMaoQYouYYangXet al. CGCG Clinical Practice Guidelines for the Management of Adult Diffuse Gliomas. Cancer Lett (2016) 375(2):263–73. doi: 10.1016/j.canlet.2016.01.024

  • CrossRef
  • Google Scholar

• 3

  LouveauASmirnovIKeyesTJEcclesJDRouhaniSJPeskeJDet al. Structural and Functional Features of Central Nervous System Lymphatic Vessels. Nature (2015) 523(7560):337–41. doi: 10.1038/nature14432

  • CrossRef
  • Google Scholar

• 4

  SahaDMartuzaRLRabkinSD. Macrophage Polarization Contributes to Glioblastoma Eradication by Combination Immunovirotherapy and Immune Checkpoint Blockade. Cancer Cell (2017) 32(2):253–67.e255. doi: 10.1016/j.ccell.2017.07.006

  • CrossRef
  • Google Scholar

• 5

  LiGWangZZhangCLiu XCaiJWangZet al. Molecular and Clinical Characterization of TIM-3 in Glioma Through 1,024 Samples. Oncoimmunology (2017) 6(8):e1328339. doi: 10.1080/2162402X.2017.1328339

  • CrossRef
  • Google Scholar

• 6

  KravchenkoJCorsiniEWilliamsMADeckerWManjiliMHOtsukiTet al. Chemical Compounds From Anthropogenic Environment and Immune Evasion Mechanisms: Potential Interactions. Carcinogenesis (2015) 36(Suppl 1):S111–127. doi: 10.1093/carcin/bgv033

  • CrossRef
  • Google Scholar

• 7

  WangZZhangCLiuXWangZSunLLiGet al. Molecular and Clinical Characterization of PD-L1 Expression at Transcriptional Level Via 976 Samples of Brain Glioma. Oncoimmunology (2016) 5(11):e1196310. doi: 10.1080/2162402X.2016.1196310

  • CrossRef
  • Google Scholar

• 8

  CaiJChenQCuiYDongJChenMWuPet al. Immune Heterogeneity and Clinicopathologic Characterization of IGFBP2 in 2447 Glioma Samples. Oncoimmunology (2018) 7(5):e1426516. doi: 10.1080/2162402X.2018.1426516

  • CrossRef
  • Google Scholar

• 9

  Eckel-PassowJELachanceDHMolinaroAMWalshKMDeckerPASicotteHet al. Glioma Groups Based on 1p/19q, IDH, and TERT Promoter Mutations in Tumors. N Engl J Med (2015) 372(26):2499–508. doi: 10.1056/NEJMoa1407279

  • CrossRef
  • Google Scholar

• 10

  VerhaakRGHoadleyKAPurdomEWangVQiYWilkersonMDet al. Integrated Genomic Analysis Identifies Clinically Relevant Subtypes of Glioblastoma Characterized by Abnormalities in PDGFRA, Idh1, EGFR, and NF1. Cancer Cell (2010) 17(1):98–110. doi: 10.1016/j.ccr.2009.12.020

  • CrossRef
  • Google Scholar

• 11

  MengXDuanCPangHChenQHanBZhaCet al. DNA Damage Repair Alterations Modulate M2 Polarization of Microglia to Remodel the Tumor Microenvironment Via the P53-Mediated MDK Expression in Glioma. EBioMedicine (2019) 41:185–99. doi: 10.1016/j.ebiom.2019.01.067

  • CrossRef
  • Google Scholar

• 12

  LimMXiaYBettegowdaCWellerM. Current State of Immunotherapy for Glioblastoma. Nat Rev Clin Oncol (2018) 15(7):422–42. doi: 10.1038/s41571-018-0003-5

  • CrossRef
  • Google Scholar

• 13

  McLendonRFriedmanABignerDFriedmanABignerDVan MeirEGet al. Comprehensive Genomic Characterization Defines Human Glioblastoma Genes and Core Pathways. Nature (2008) 455(7216):1061–8. doi: 10.1038/nature07385

  • CrossRef
  • Google Scholar

• 14

  WiendlHMitsdoerfferMHofmeisterVWischhusenJBornemannAMeyermannRet al. A Functional Role of HLA-G Expression in Human Gliomas: An Alternative Strategy of Immune Escape. J Immunol (2002) 168(9):4772–80. doi: 10.4049/jimmunol.168.9.4772

  • CrossRef
  • Google Scholar

• 15

  DidenkoVVNgoHNMinchewCBaskinDS. Apoptosis of T Lymphocytes Invading Glioblastomas Multiforme: A Possible Tumor Defense Mechanism. J Neurosurg (2002) 96(3):580–4. doi: 10.3171/jns.2002.96.3.0580

  • CrossRef
  • Google Scholar

• 16

  ParsaATWaldronJSPannerACraneCAParneyIFBarryJJet al. Loss of Tumor Suppressor PTEN Function Increases B7-H1 Expression and Immunoresistance in Glioma. Nat Med (2007) 13(1):84–8. doi: 10.1038/nm1517

  • CrossRef
  • Google Scholar

• 17

  WherryEJ. T Cell Exhaustion. Nat Immunol (2011) 12(6):492–9. doi: 10.1038/ni.2035

  • CrossRef
  • Google Scholar

• 18

  GreterMHeppnerFLLemosMPOdermattBMGoebelsNLauferTet al. Dendritic Cells Permit Immune Invasion of the CNS in an Animal Model of Multiple Sclerosis. Nat Med (2005) 11(3):328–34. doi: 10.1038/nm1197

  • CrossRef
  • Google Scholar

• 19

  YoshiharaKShahmoradgoliMMartinezEVegesnaRKimHTorres-GarciaWet al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

  • CrossRef
  • Google Scholar

• 20

  VerhaakRGTamayoPYangJYHubbardDZhangHCreightonCJet al. Prognostically Relevant Gene Signatures of High-Grade Serous Ovarian Carcinoma. J Clin Invest (2013) 123(1):517–25. doi: 10.1172/JCI65833

  • CrossRef
  • Google Scholar

• 21

  YangJZhangDFrangiAFYangJY. Two-Dimensional PCA: A New Approach to Appearance-Based Face Representation and Recognition. IEEE Trans Pattern Anal Mach intell (2004) 26(1):131–7. doi: 10.1109/TPAMI.2004.1261097

  • CrossRef
  • Google Scholar

• 22

  MayakondaALinDCAssenovYPlassCKoefflerHP. Maftools: Efficient and Comprehensive Analysis of Somatic Variants in Cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118

  • CrossRef
  • Google Scholar

• 23

  NewmanAMLiuCLGreenMRGentlesAJFengWXuYet al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

  • CrossRef
  • Google Scholar

• 24

  ThommenDSKoelzerVHHerzigPRollerATrefnyMDimeloeSet al. A Transcriptionally and Functionally Distinct PD-1(+) Cd8(+) T Cell Pool With Predictive Potential in non-Small-Cell Lung Cancer Treated With PD-1 Blockade. Nat Med (2018) 24(7):994–1004. doi: 10.1038/s41591-018-0057-z

  • CrossRef
  • Google Scholar

• 25

  CaiMCZhaoXCaoMMaPChenMWuJet al. T-Cell Exhaustion Interrelates With Immune Cytolytic Activity to Shape the Inflamed Tumor Microenvironment. J pathol (2020) 251(2):147–59. doi: 10.1002/path.5435

  • CrossRef
  • Google Scholar

• 26

  SubramanianATamayoPMoothaVKMukherjeeSEbertBLGilletteMAet al. Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

  • CrossRef
  • Google Scholar

• 27

  HänzelmannSCasteloRGuinneyJ. Gsva: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinf (2013) 14(1):7. doi: 10.1186/1471-2105-14-7

  • CrossRef
  • Google Scholar

• 28

  GuZLeiGRolandEMatthiasSBenediktB. N Circlize\N Implements and Enhances Circular Visualization in R\N. Bioinformatics (2014) 19:19. doi: 10.1093/bioinformatics/btu393

  • CrossRef
  • Google Scholar

• 29

  KrzywinskiMScheinJBirolIConnorsJGascoyneRHorsmanDet al. Circos: An Information Aesthetic for Comparative Genomics. Genome Res (2009) 19(9):1639–45. doi: 10.1101/gr.092759.109

  • CrossRef
  • Google Scholar

• 30

  RiehmannPHanflerMFroehlichB. Interactive Sankey Diagrams. Paper Presented at: IEEE Symposium on Information Visualization, 2005. INFOVIS 2005 (2005) 2005.

  • Google Scholar

• 31

  YangWSoaresJGreningerPEdelmanEJLightfootHForbesSet al. Genomics of Drug Sensitivity in Cancer (GDSC): A Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res (2013) 41(Database issue):D955–961. doi: 10.1093/nar/gks1111

  • CrossRef
  • Google Scholar

• 32

  GeeleherPCoxNJHuangRS. Clinical Drug Response can be Predicted Using Baseline Gene Expression Levels and In Vitro Drug Sensitivity in Cell Lines. Genome Biol (2014) 15(3):R47. doi: 10.1186/gb-2014-15-3-r47

  • CrossRef
  • Google Scholar

• 33

  IorioFKnijnenburgTAVisDJBignellGRMendenMPSchubertMet al. A Landscape of Pharmacogenomic Interactions in Cancer. Cell (2016) 166(3):740–54. doi: 10.1016/j.cell.2016.06.017

  • CrossRef
  • Google Scholar

• 34

  PopeWBSayreJPerlinaAVillablancaJPMischelPSCloughesyTF. MR Imaging Correlates of Survival in Patients With High-Grade Gliomas. AJNR Am J neuroradiol (2005) 26(10):2466–74.

  • Google Scholar

• 35

  BruzzoneMGEoliMCuccariniVGrisoliMVallettaLFinocchiaroG. Genetic Signature of Adult Gliomas and Correlation With MRI Features. Expert Rev Mol diagn (2009) 9(7):709–20. doi: 10.1586/erm.09.44

  • CrossRef
  • Google Scholar

• 36

  CarrilloJALaiANghiemphuPLPhillipsHSKharbandaSMoftakharPet al. Relationship Between Tumor Enhancement, Edema, IDH1 Mutational Status, MGMT Promoter Methylation, and Survival in Glioblastoma. AJNR Am J neuroradiol (2012) 33(7):1349–55. doi: 10.3174/ajnr.A2950

  • CrossRef
  • Google Scholar

• 37

  ThorssonVGibbsDLBrownSDWolfDBortoneDSYangT-HOet al. The Immune Landscape of Cancer. Immunity (2018) 48(4):812–30.e814. doi: 10.1016/j.immuni.2018.03.023

  • CrossRef
  • Google Scholar

• 38

  LiuZLiMJiangZWangX. A Comprehensive Immunologic Portrait of Triple-Negative Breast Cancer. Trans Oncol (2018) 11(2):311–29. doi: 10.1016/j.tranon.2018.01.011

  • CrossRef
  • Google Scholar

• 39

  JiangZLiuZLiMChenCWangX. Immunogenomics Analysis Reveals That TP53 Mutations Inhibit Tumor Immunity in Gastric Cancer. Trans Oncol (2018) 11(5):1171–87. doi: 10.1016/j.tranon.2018.07.012

  • CrossRef
  • Google Scholar

• 40

  LyuHLiMJiangZLiuZWangX. Correlate the TP53 Mutation and the HRAS Mutation With Immune Signatures in Head and Neck Squamous Cell Cancer. Comput Struct Biotechnol J (2019) 17:1020–30. doi: 10.1016/j.csbj.2019.07.009

  • CrossRef
  • Google Scholar

• 41

  NakashimaHAlayoQAPenaloza-MacMasterPFreemanGJKuchrooVKReardonDAet al. Modeling Tumor Immunity of Mouse Glioblastoma by Exhausted CD8(+) T Cells. Sci Rep (2018) 8(1):208–8. doi: 10.1038/s41598-017-18540-2

  • CrossRef
  • Google Scholar

• 42

  Pombo AntunesARScheyltjensIDuerinckJNeynsBMovahediKVan GinderachterJA. Understanding the Glioblastoma Immune Microenvironment as Basis for the Development of New Immunotherapeutic Strategies. eLife (2020) 9. doi: 10.7554/eLife.52176

  • CrossRef
  • Google Scholar

• 43

  MummJBEmmerichJZhangXChanIWuLMauzeSet al. Il-10 Elicits IFNgamma-dependent Tumor Immune Surveillance. Cancer Cell (2011) 20(6):781–96. doi: 10.1016/j.ccr.2011.11.003

  • CrossRef
  • Google Scholar

• 44

  MantovaniAMarchesiFMalesciALaghiLAllavenaP. Tumour-Associated Macrophages as Treatment Targets in Oncology. Nat Rev Clin Oncol (2017) 14(7):399–416. doi: 10.1038/nrclinonc.2016.217

  • CrossRef
  • Google Scholar

• 45

  KlemmFMaasRRBowmanRLKorneteMSoukupKNassiriSet al. Interrogation of the Microenvironmental Landscape in Brain Tumors Reveals Disease-Specific Alterations of Immune Cells. Cell (2020) 181(7):1643–60.e1617. doi: 10.1016/j.cell.2020.05.007

  • CrossRef
  • Google Scholar

• 46

  HugoWZaretskyJMSunLSongCMorenoBHHu-LieskovanSet al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell (2016) 165(1):35–44. doi: 10.1016/j.cell.2016.02.065

  • CrossRef
  • Google Scholar

• 47

  ZhaCMengXLiLMiSQianDLiZet al. Neutrophil Extracellular Traps Mediate the Crosstalk Between Glioma Progression and the Tumor Microenvironment Via the HMGB1/RAGE/IL-8 Axis. Cancer Biol Med (2020) 17(1):154–68. doi: 10.20892/j.issn.2095-3941.2019.0353

  • CrossRef
  • Google Scholar

• 48

  RodriguesJCGonzalezGCZhangLIbrahimGKellyJJGustafsonMPet al. Normal Human Monocytes Exposed to Glioma Cells Acquire Myeloid-Derived Suppressor Cell-Like Properties. Neuro-Oncology (2010) 12(4):351–65. doi: 10.1093/neuonc/nop023

  • CrossRef
  • Google Scholar

• 49

  CipollettaDFeuereMLiAKameiNLeeJShoelsonSEet al. Ppar-γ is a Major Driver of the Accumulation and Phenotype of Adipose Tissue T_(reg) Cells. Nature (2012) 486(7404):549–53. doi: 10.1038/nature11132

  • CrossRef
  • Google Scholar

• 50

  RangachariMZhuCSakuishiKXiaoSKarmanJChenAet al. Bat3 Promotes T Cell Responses and Autoimmunity by Repressing Tim-3-mediated Cell Death and Exhaustion. Nat Med (2012) 18(9):1394–400. doi: 10.1038/nm.2871

  • CrossRef
  • Google Scholar

• 51

  ElTanboulyMAZhaoYNowakELiJSchaafsmaEMercierILet al. VISTA is a Checkpoint Regulator for Naïve T Cell Quiescence and Peripheral Tolerance. Science (2020) 367(6475):eaay0524. doi: 10.1126/science.aay0524

  • CrossRef
  • Google Scholar

• 52

  TopalianSLTaubeJMPardollDM. Neoadjuvant Checkpoint Blockade for Cancer Immunotherapy. Science (2020) 367(6477):eaax0182. doi: 10.1126/science.aax0182

  • CrossRef
  • Google Scholar

• 53

  ImSJHashimotoMGernerMYLeeJKissickHTBurgerMCet al. Defining CD8+ T Cells That Provide the Proliferative Burst After PD-1 Therapy. Nature (2016) 537(7620):417–21. doi: 10.1038/nature19330

  • CrossRef
  • Google Scholar

• 54

  LiBSeversonEPignonJCZhaoHLiTNovakJet al. Comprehensive Analyses of Tumor Immunity: Implications for Cancer Immunotherapy. Genome Biol (2016) 17(1):174. doi: 10.1186/s13059-016-1028-7

  • CrossRef
  • Google Scholar

• 55

  ChenQHanBMengXDuanCYangCWuZet al. Immunogenomic Analysis Reveals LGALS1 Contributes to the Immune Heterogeneity and Immunosuppression in Glioma. Int J Cancer (2019) 145(2):517–30. doi: 10.1002/ijc.32102

  • CrossRef
  • Google Scholar

• 56

  SanmamedMFChenL. A Paradigm Shift in Cancer Immunotherapy: From Enhancement to Normalization. Cell (2018) 175(2):313–26. doi: 10.1016/j.cell.2018.09.035

  • CrossRef
  • Google Scholar

• 57

  JinXKimLJYWuQWallaceLCPragerBCSanvoranartTet al. Targeting Glioma Stem Cells Through Combined BMI1 and EZH2 Inhibition. Nat Med (2017) 23(11):1352–61. doi: 10.1038/nm.4415

  • CrossRef
  • Google Scholar

• 58

  YarchoanMHopkinsAJaffeeEM. Tumor Mutational Burden and Response Rate to PD-1 Inhibition. N Engl J Med (2017) 377(25):2500–1. doi: 10.1056/NEJMc1713444

  • CrossRef
  • Google Scholar

• 59

  CampbellPJGetzGKorbelJOStuartJMJenningsJLSteinLDet al. Pan-Cancer Analysis of Whole Genomes. Nature (2020) 578(7793):82–93. doi: 10.1038/s41586-020-1969-6

  • CrossRef
  • Google Scholar