The genetically altered landscape of immune-evasion regulators in LGGs
The workflow of this study is depicted in Fig. 1. To assess the genomic alterations of IEVGs in LGGs, the transcriptional expression, copy number variations (CNVs) and single nucleotide polymorphisms (SNPs) of 179 cancer-intrinsic immune evasion genes were evaluated (3 of 182 IEVGs were unmapped in the TCGA cohort). Enrichment analysis of these 179 IEVGs using the Metascape webtool showed that, when GO-biological processes and KEGG pathway terms were included, these IEVGs were mainly enriched in the tumor necrosis factor (TNF) signaling pathway, autophagy, ferroptosis, antigen processing and presentation, ubiquitin mediated proteolysis, Ras signaling pathway, apoptosis, and protein processing in the endoplasmic reticulum (Fig. 2A, B). Several of these processes and pathways were previously shown to be essential in regulating cancer immune evasion [41,42,43]. The overall differences of transcriptional expression between LGGs and normal brain tissues (NBTs), including the CNV and SNP proportions of each IEVG in LGGs, were visualized using a circus heatmap (Fig. 2C). In this circus heatmap, the location of each IEVG is depicted on chromosomes, and the expression levels compared with NBTs, CNVs (including diploid, gain and loss) and SNP proportions are shown in the inner heatmaps. Blue indicates high expression or proportion, whereas red indicates lower expression or proportion. Most of the 179 IEVGs were upregulated in LGGs compared with NBTs, with the top 20 upregulated IEVGs shown in the boxplots (Fig. 2D). Ninety-three IEVGs were significantly overexpressed in LGGs, with the most upregulated IEVG being CEP55. The top 20 IEVGs with most frequent CNVs are shown in Fig. 2E. The alterations in CNV proportions of 14 IEVGs were > 30%, whereas the alterations of 130 IEVGs were > 10%. GPI was the IEVG with most frequent CNVs in LGGs (47.56%; gain 41.52%; loss, 6.04%). Analysis of mutation frequency of IEVGs in LGGs showed that all of the IEVGs had mutation frequencies < 1% (Fig. 2F), indicating that mutations of IEVGs had little effect on LGG immune evasion. Taken together, these findings show that aberrant genomic alterations were present in 172 IEVGs of LGGs, with CNV being the primary type of IEVG alterations. Because CEP55 was the most significantly upregulated gene in LGGs, the expression of CEP55 protein was analyzed in seven paired LGG cores and adjacent tissues (Fig. 2G, full-length gels are presented in Figure S1), with these findings showing that the CEP55 protein was also upregulated in LGGs (Fig. 2H).
Identification of intrinsic immune-evasion patterns in LGGs
The cancer-intrinsic immune-evasion patterns in the LGG TME were assessed by unsupervised clustering of LGG patients based on the transcriptional profile of 88 prognostic IEVGs that had been identified by univariate Cox regression analysis. Clustering analysis showed that k = 3 was the optimal cluster number in LGG samples (Figure S2A-D). Clustering analysis divided LGG patients in the TCGA-LGG cohort into three distinct clusters, with Clusters 1, 2, and 3 consisting of 175, 104, and 178 patients, respectively. Unsupervised clustering of the two merged LGG cohorts, using k = 3 as the optimal clustering number k = 3 was also performed. In order to investigate the potential clinical implications of the IEVG-based cluster, overall survival (OS) was compared in each pair of clusters in each LGG cohort. In the TCGA cohort, LGG patients in Cluster 2 had the poorest clinical prognosis when compared with patients in Clusters 1 (p = 0.0002) and 3 (p < 0.0001), whereas the difference between Clusters 1 and 3 was not statistically significant (p = 0.076) (Fig. 3A). Kaplan–Meier survival analysis in the meta-GEO yielded consistent results, indicating that OS was poorest in Cluster 2 (Fig. 3B). Further analysis of the meta-CGGA cohort showed that patients in Cluster 3 had significantly better prognosis than patients in Clusters 1 (p = 0.00028) and 2 (p < 0.0001), with the difference between Clusters 1 and 2 also being significant (p < 0.0001) (Fig. 3C). These findings indicated that the clusters might lead to promising clinical applications in patients with LGG.
To further analyze the underlying differences in biological processes and/or pathways between Cluster 2 and the other two clusters, GSVA scoring of hallmarks for LGG patients were compared in these clusters. Immune-activated associated processes, including interferon-γ, interferon-α, and inflammatory responses and TNFA signaling via NFKB, were significantly enriched in Cluster 2 compared with Clusters 1 and 3 (Fig. 3D, E). In addition, previously described signature scores were compared in the three IE clusters (Fig. 3F), with signature scores for antigen processing, CD8 T effector and immune checkpoint being much higher in Cluster 2 than in Clusters 1 and 3, similar to the GSVA results. Other characteristics of malignancy, including angiogenesis, cell cycle associated processes, EMT markers and DNA repair signatures, were also significantly higher in Cluster 2, whereas stromal associated signatures, including FGFR3 related genes and WNT target, were lowest in Cluster 2. These findings showed that the application of unsupervised clustering based on the profiles of IEVGs in LGGs could classify LGG into three subgroups, with distinct clinical prognosis and tumor immune phenotypic signatures. The abundance of immune-cell infiltration was subsequently incorporated to show infiltration differences among the three clusters.
Associations between LGG cancer-intrinsic immune-evasion patterns and distinct immune cell infiltration
To further evaluate the molecular characteristics and immune cell infiltration among the three clusters, IDH mutations and 1p/19q co-deletion status were compared. Significant differences were observed among the three clusters, with most LGGs in Cluster 1 containing IDH mutations and 1p/19q non-codeletions, most LGGs in Cluster 2 containing IDH-wild type, and most LGGs in Cluster 3 containing IDH-mutations and 1p/19q codeletions (Fig. 4A). LGGs in Cluster 2 had the highest ESTIMATE score, which correlated positively with immune cell infiltration level, and the lowest tumor purity, whereas LGGs in Cluster 3 had the lowest ESTIMATE score and the highest tumor purity (Fig. 4B, C). These findings indicated that LGGs in Cluster 2 were infiltrated by larger numbers of immune and/or stromal cells, whereas LGGs in Cluster 3 LGGs were not infiltrated by extraneous cells, with LGGs in Cluster 1 showing some immune cell infiltration.
Principal component analysis (PCA) showed that the three LGG Clusters occupied distinct spaces in the three-dimension PCA model based on whole transcriptional profiles (Fig. 4D), indicating that transcriptional discrepancies may lead to differences in immune-evasion among LGGs. Quantification of the abundance of 28 immune cell types in each LGG sample using the ssGSEA algorithm showed significant differences in immune cell infiltrations among the three LGG Clusters (Fig. 4E). CIBERSORT quantification of the levels of infiltration of various immune cells showed that CD8 + T cells, CD4 + memory resting T cells, regulatory T cells and M2 macrophages were enriched in the LGGs in Cluster 2, whereas follicular helper T cells were enriched in Cluster 3 (Fig. 4F).
Identifying DEGs among the three LGG clusters
Although LGGs were clustered into three immune-evasion phenotypes, the potential differences in transcriptional alteration among these clusters were not clearly understood. The differences in transcriptional expression between each pair of clusters were investigated using a Bayesian method to recognize overlapping differentially expressed genes (DEGs). A total of 607 DEGs that represented crucial distinguishing scores of the three LGG clusters were considered IEV signature genes (Fig. 5A). Unsupervised consensus cluster analysis of the 607 DEGs in the TCGA-LGG cohort yielded three novel stable transcriptomic LGG clusters (Figure S3 A–D; k = 3), with a heatmap showing the expression patterns of the 607 DEGs in the three clusters (Fig. 5B). Gene ontology biological processes (GO-BP) analysis of the 607 DEGs showed enrichment of many onco-immunology related processes, including T cell activation, response to interferon-γ, neutrophil mediated immunity, and regulation of immune effector processes (Fig. 5C). These results further indicated that the 607 DEGs were associated with immune evasion processes and that the novel classifications might have clinical implications. Evaluation of survival showed that LGGs in Cluster S2 were associated poorer prognosis than the LGGs in the other two clusters (Fig. 5D). Moreover, univariate Cox regression analysis indicated that Cluster S2 was associated with a significant risk than Cluster S1 (HR = 3.28; 95% CI: 2.00–5.36, p < 0.001), whereas risks did not differ significantly in Clusters S3 and Cluster S1 (HR = 0.64; 95% CI: 0.34–1.20, p = 0.164) (Fig. 5E). Comparisons of the levels of 19 biological processes and pathway signatures showed that LGGs in Cluster S2 were significantly enriched in immune response processes, including antigen processing, CD8 T effector, immune checkpoint, EMT associated signatures (EMT1/2/3), cell cycle and DNA repair related signatures, including cell cycle, DNA damage repair, DNA replication, homologous recombination and mismatch repair, whereas FGFR3 related genes and WNT target signatures were enriched in Clusters 1 and 3.
Determination of IEVscores and their clinical implications
Although these findings showed that immune-evasion patterns in LGGs had prognostic implications and were associated with immune cell infiltration, a method of evaluating the immune-evasion level in LGG patients has not yet been determined. IEVscore, a method based on the levels of expression of IE signature genes, was used to quantify the degree of immune evasion in individual patients. The transcriptomic expression profiles of the 607 DEGs were subjected to univariate Cox regression analysis, followed by Lasso regression analysis. These processes identified 17 IESig genes that were used to establish an IEVscore scoring system (Figure S4A, B). The workflow involved in constructing the IEVscore is shown in a Sankey plot (Fig. 6A). To better understand the underlying biological processes and pathways associated with the IEVscore, the statistical correlations between IEVscore and the 19 signatures were determined by Spearman correlation analysis. The resulting correlation heatmap revealed that the IEVscore was positively correlated with the signatures for antigen processing, CD8 T effector and immune checkpoint (Fig. 6B), indicating that higher IEVscores might indicate a higher degree of immune activation in the microenvironment of LGGs.
The ROC curves revealed that IEVscore was a highly robust predictor of OS in LGG patients, with the AUC values of IEVscores predicting 1–5-year OS were > 0.6 in all three patient cohorts (Fig. 6C). Survival analysis indicated that higher IEV scores in all three LGG cohorts were significantly prognostic of poorer outcomes (Fig. 6D–F, p < 0.0001). The associations between IEVscores and the molecular subtypes of LGGs were assessed by comparison analysis, with the results showing that IDH-wild type LGGs were associated with significantly higher IEVscores in all three cohorts, with IDH-mut and 1p/19q co-deletion LGGs having the lowest IEVscores in the TCGA and meta-CGGA cohorts (Fig. 6G-I). Because IEVscore may play a role in response to immunotherapy, the tumor mutation burden (TMB) level was compared in the low- and high-IEVscore subgroups in the TCGA cohort. These comparisons showed that TMB level was significantly higher in the high- than in the low-IEVscore subgroup (Fig. 6J), indicating that IEVscore not only predicts the prognosis of patients with LGG but may also reflect responses to immunotherapy. Furthermore, a waterfall plot comparing the mutational landscapes in the low- and high-IEVscore subgroups in the TCGA-cohort showed the most frequently mutated genes in these subgroups (Fig. 6K, L).
IEVscore predicts response to immunotherapy
To determine whether IEVscore could predict response to cancer immunotherapy, TIDE scores were calculated, and the responses of LGG patients to immunotherapy were predicted using the TIDE algorithm. LGGs with higher IEVscores also showed higher TIDE scores and fewer predicted responders to immunotherapy in all the three LGG cohorts (Fig. 7A–C). Further, IEVscores not only predicted the prognosis of cancer patients treated with immune checkpoint inhibitors (ICIs) but also reflected the immunotherapy response status of cancer patients. Kaplan–Meier survival analysis showed that IEVscore stratified melanoma patients (Vanallen cohort) into two groups with different prognoses, with patients having high-IEVscores showing poorer clinical outcomes (Fig. 7D) and greater resistance to anti-CTLA4 therapy (Fig. 7E). Cancer patients in the IMvigor210 cohort treated with anti-PDL1 could be classified into high- and low-IEVscore groups, with high-IEV score being associated with shorter OS (Fig. 7F) and a lower percentage responding to ICIs (Fig. 7G).
TMZ sensitivity analysis
Because TMZ is the primary postoperative chemotherapy agent currently used to treat patients with LGG, the correlation between TMZ sensitivity and IEVscore was investigated in these patients. The relative AUC values representing the levels of TMZ resistance were obtained by ridge-regression for the 1421 LGG patients in the Cancer Therapeutics Response Portal (CTRP) and PRISM Repurposing datasets. Patients with high-IEVscore LGG showed greater TMZ resistance in all three LGG cohorts using both the PRISM and CTRP algorithms (Fig. 7H–J). These findings indicated that LGG patients with lower IEVscore might benefit more from TMZ therapy. Subsequent, the prognostic predictive power of IEVscore was investigated in LGG patients in the TCGA and meta-CGGA datasets who were and were not treated with TMZ. Kaplan–Meier survival analysis indicated that IEVscore remained prognostic role and stratified patients who were or were not treated with TMZ into subgroups with distinct prognosis (Fig. 7K–M). Furthermore, Pearson correlation analysis showed that the AUC values, which represent the drug resistance levels of cancer patients, were significantly positively correlated with IEscores in all three LGG cohorts by using the CTRP (Figure S5 A–C) and PRISM (Figure S5 D–F) databases.
