
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)13017-4
10.1016/j.heliyon.2024.e36986
e36986
Research Article
Molecular subtypes based on immunologic and epithelial-mesenchymal transition gene sets reveal tumor immune microenvironment characteristics and implications for immunotherapy of patients with glioma
Lin Wen-wen 15wwlin@alumni.stu.edu.cn
a⁎1
Zhao Wei-jiang b1
Ou Guan-yong 12231350@mail.sustech.edu.cn
c⁎⁎
a Department of Pathology, Shenzhen Hospital, Southern Medical University, Shenzhen, 518101, China
b Cell biology Department, Wuxi School of Medicine, Jiangnan University, Wuxi, 214122, China
c School of Medicine, Southern University of Science and Technology, Shenzhen, 518055, China
⁎ Corresponding author. 15wwlin@alumni.stu.edu.cn
⁎⁎ Corresponding author. 12231350@mail.sustech.edu.cn
1 Contributed equally.

27 8 2024
15 9 2024
27 8 2024
10 17 e3698628 11 2023
24 8 2024
26 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
The tumor immune microenvironment (TIME) significantly influences cancer progression and treatment. This study sought to uncover novel TIME-related glioma biomarkers to advance antitumor immunotherapies by integrating data from sequencing of bulk RNA as well as scRNA. Immunologic and epithelial-mesenchymal transition (EMT) characteristics were used to classify glioma patients into two immune subtypes (ISs) and two EMT subtypes (ESs). Patients in IS1 and ES1, characterized by high immune infiltration and low stemness scores, exhibited poor clinical outcomes and limited responsiveness to immunotherapy. A new risk signature was developed using 16 genes and validated in independent glioma cohorts. Among these, HAVCR2, IL18, LAGLS9, and PTPN6 emerged as hub genes, with IL18 identified as a potential independent indicator. The upregulation of IL18 in high-grade gliomas and U-251 MG cells aligned with bioinformatics analysis. These insights deepen the understanding of TIME-related mechanisms in glioma and highlight potential therapeutic targets, offering a theoretical foundation for effective antitumor immunotherapies in glioma.

Keywords

Molecular subtypes
Immunologic and epithelial-mesenchymal transition gene sets
Tumor immune microenvironment
Immunotherapy
==== Body
pmc1 Introduction

Gliomas exist as the most common brain tumor that arises from cells in the glia, significantly contributing to morbidity and mortality [1]. Glioblastoma (GBM) stands as the most frequent variant of malignant glioma, comprising roughly 80 % of primary malignant brain tumors. It has a grim prognosis, with under two years of median overall survival (OS), significantly diminishing the quality of life [2]. Despite advances in understanding glioma biology and the development of targeted therapies in preclinical and clinical trials, the incidence of recurrence, cognitive deficiencies, disability, and mortality associated with glioma remains high [3]. Most glioma treatments with notable survival benefits, such as radiation and chemotherapy, rely on nonspecific targeting of proliferating cells and are often compromised by glioma cell invasion and immune evasion [3,4]. Therefore, it is crucial to identify unknown factors or further investigate known factors that influence the tumor microenvironment (TME) to promote glioma progression, malignant transformation, and drug resistance. This approach is vital for advancing the assessment and customization of therapeutic options for gliomas.

The TME, which is the site of tumor emergence, is significantly implicated in its development [5,6]. TME primarily consists of the extracellular matrix, various cellular components such as fibroblasts and immune components, along with vascular networks that are sourced from adjoining tissues [[5], [6], [7]]. With the participation of immune-related genes (IRGs) and epithelial-mesenchymal transition (EMT)-related genes (ERGs), the TME functions crucially in tumor progression and treatment resistance, influenced by extracellular matrices, cytokines, growth signals, and other elements, thereby affecting cancer treatment [8].

Compared to brain tissue in the normal state, which exhibits a low to moderate immune response, the tumor immune microenvironment (TIME) significantly contributes to the aggressive advancement and therapeutic resistance in glioma [9]. A suppressive TIME can severely impair GBM treatment [10]. Remodeling TIME using a liposomal honokiol and disulfiram/copper codelivery system (CDX-LIPO) specifically targeting the brain can trigger autophagy and immune-mediated cell death in tumor cells. This process activates macrophages and dendritic cells infiltrating tumors, as well as the priming of T and NK (natural killer) cells, leading to an immunity response against tumors and tumor shrinkage [10]. Recently, various IRGs have been identified, which can be categorized as low-risk and high-risk [[11], [12], [13]], suggesting that IRGs play differential and complex roles in the emergence, advancement, and metastatic spread of gliomas.

Increasing evidence suggests EMT facilitates glioma invasion and progression. EMT, a multifactor-regulated process, involves various induction agents, families of transcription factors, and an extensive range of genes involved in signaling pathways [14]. Enhanced EMT is associated with poor survival, whereas the EMT in the initial stages iscomparatively less malignant in GBM [15], suggesting that EMT assessment could predict GBM prognosis and aid in the advancement of innovative therapeutic strategies and prognostic markers [15]. Additionally, the elevated expression of EMT genes in GBMs and pilocytic astrocytomas is indicative of mural cells proliferation within newly formed blood vessels [16].

Six collagen genes have been identified as regulators of both the glioma immunosuppressive tumor microenvironment and the EMT process, namely COL1A1, COL1A2, COL3A1, COL4A1, COL4A2, and COL5A2 [17]. This suggests that IRGs and ERGs may synergistically influence glioma behaviors such as invasion, migration, and metastasis. Cohort analysis of The Cancer Genome Atlas (TCGA) uncovered two immune subtypes (ISs) and two EMT subtypes (ESs) in patients with glioma based on IRG and ERG sets. This resulted in the elucidation of key immune- and EMT-related genes (IERGs), with sixteen chosen to construct a novel risk prediction model for glioma patients. Among these, four IERGs (HAVCR2, IL18, LAGLS9, and PTPN6) emerged as hub genes, showing high expression in M1 macrophages and monocytes. Furthermore, the upregulation of IL18 in GBM indicates its potential as an independent clinical marker for patients with high-grade glioma. This integrated analysis of IRGs and ERGs in gliomas establishes a theoretical framework aimed at developing and implementing more effective antitumor immunotherapies.

2 Results

2.1 Identification of novel ISs and ESs of gliomas and analysis of their survival and different immune features

The study framework is detailed in Supplementary Fig. 1. Consensus clustering without supervision (k = 2) was applied to 663 gliomas from the TCGA cohort, focusing on immunologic and EMT gene signatures. This analysis categorized patients with glioma into two ISs (IS1 and IS2) and two ESs (ES1 and ES2) (Fig. 1A–H). Based on Kaplan-Meier survival analysis, individuals in the IS1 and ES1 subtypes had better OS compared to those in IS2 and ES2 (Fig. 1B–I). Gene set enrichment analysis (GESA) of HALLMARK pathways revealed significant differential enrichment in five pathways—interferon-gamma response, IL2-STAT5 signaling, apoptosis, hypoxia, and glycolysis—in IS1 (Fig. 1C) and ES1 (Fig. 1J) based on NES, P-value, and FDR value.Fig. 1 Unsupervised consensus clustering of gliomas based on immunologic and EMT gene signatures in the TCGA cohort. (A and H) Heatmap depicting consensus clustering matrix (k = 2) for 1959 immunologic genes (A) and 1263 EMT genes (H) across 663 glioma samples. (B and I) Kaplan-Meier survival analysis of patients with gliomas exhibiting different ISs (B) and ESs (I). (C and J) GSEA for pathways in HALLMARK terms of different ISs (C) and ESs (J). (D and K) Potential ICB response predicted with TIDE algorithm in different ISs (D) and ESs (K). (E and L) Stemness score based on OCLR algorithm in different ISs (E) and ESs (L). (F-G and M-N) TMB and MSI score in different ISs (F-G) and ESs (M-N).

Fig. 1

The TIDE algorithm, used to predict potential immunotherapeutic responses, showed significantly higher TIDE scores in IS1 and ES1 compared to IS2 and ES2 (Fig. 1D–K), indicating that IS1 and ES1 patients had poorer responses to immune checkpoint blockade (ICB) immunotherapy and experienced decreased survival following therapies. The OCLR algorithm assessed stemness scores, showing higher scores in IS2 and ES2 than in IS1 and ES1 (Fig. 1E–L), which correlates with greater tumor differentiation in IS2 and ES2 patients. Tumor mutational burden (TMB) and microsatellite instability (MSI) scores, critical for predicting immune efficacy, were higher in IS1 and ES1 compared to IS2 and ES2 (Fig. 1F–M) and lower in MSI score (Fig. 1G–N). Interestingly, IS1 and ES1 exhibited similar survival prognoses and clinicopathological features.

2.2 TIME features of distinct TCGA cohort subtypes

To evaluate the association of glioma subtypes within immune scores, the ESTIMATE algorithm evaluated stromal and immune cell infiltration, along with tumor purity, across different ISs and ESs. This analysis generated four scores: immune, stromal, ESTIMATE, and tumor purity. The results indicated IS1 and ES1 showed increased immune scores (Supplemental Figs. 2A and E), stromal scores (Supplemental Figs. 2B and F), and ESTIMATE scores (Supplemental Figs. 2C and G), along with lower tumor purity (Supplemental Figs. 2D and H) compared to IS2 and ES2. Higher tumor purity in IS2 and ES2 suggested fewer infiltrated immune and stromal cells in these gliomas. The infiltration profiles of 22 distinct immune cell types were assessed by the CIBERSORT algorithm across various ISs and ESs, and the results were displayed using heatmaps (Supplemental Figs. 2I and K). Clear differences were noted in various populations of immune cells in the comparison of IS1 and IS2, and between ES1 and ES2. Notably, the composition of B cell plasma, resting T cells CD4+ memory, monocytes, and M2 macrophages varied significantly between IS1 and IS2 and between ES1 and ES2 (Supplemental Figs. 2J and L).

2.3 Linkage between ISs and ESs with ICPs, ICD modulators, and HLA family

Considering the critical roles of immune checkpoint proteins (ICPs), immunogenic cell death (ICD) modulators, and HLA family genes in cancer immunity and the effectiveness of vaccines against mRNA, distinct expression of these elements was assessed in different ISs and ESs. The levels of expression for 25 ICD modulators, 47 ICPs, and 24 HLA family genes were investigated across different ISs and ESs. The analysis revealed that 19 ICD modulators, 39 ICPs, and 23 HLA-related genes were differentially expressed in these subtypes (Supplemental Fig. 3). For instance, 18 ICD modulators were significantly upregulated in IS1 and ES1 gliomas. Additionally, 38 ICPs, excluding ADORA2A, CD27, IDO2, KIR3DL1, TIGIT, TNFRSF25, TNFSF18, TNFSF9, and VSIR, showed significant upregulation in IS1 and ES1 gliomas (Supplemental Fig. 3). All HLA-related genes, except HLA.J, exhibited significantly higher expression levels in IS1 and ES1 gliomas (Supplemental Fig. 3). These findings suggest that IS1 and ES1 gliomas share similar expression profiles of ICPs, ICD modulators, and HLA family genes, indicating that patients with IS1 and ES1 gliomas may benefit from similar mRNA vaccine therapies.

2.4 Exploring Time-related modules by WGCNA

To recognize IERGs, WGCNA on the gliomas-TCGA cohort constructed co-expression networks and pinpointed gene modules linked with immune scores. Setting the soft threshold to 12 (Fig. 2A), the analysis identified 20 key TIME-related gene modules, each color-coded (Fig. 2B). The eigengene adjacency heatmap illustrated correlations among these modules (Fig. 2C). Further exploration revealed their associations with clinical factors, such as grade, age, gender, and various immune-related parameters (Fig. 2D). The brown module, showing the strongest links to stromal, immune and ESTIMATE score, was selected to conduct deeper analysis (Fig. 2D). Scatter plots demonstrated significant positive associations of the brown module membership with immune score (Fig. 2E), stromal score (Fig. 2F), ESTIMATE score (Fig. 2G), and tumor purity (Fig. 2H).Fig. 2 WGCNA screening for TIME-related module in the TCGA cohort. (A) Analysis of the scale-free topology fit index and mean connectivity for various soft-threshold powers (β), selecting a power of 12. (B) Dendrogram of all genes clustered based on a dissimilarity measure (1-TOM). (C) Correlation heatmap of different modules. (D) Determination of module-trait relationships in patients with glioma. (E-H) Scatter plots of module eigengene significance for immune score (E), stromal score (F), ESTIMATE score (G), and tumor purity (H) in the brown module. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 2

2.5 Identification of IERGs and functional enrichment analysis

The upset plot identified four hub intersection genes (HAVCR2, IL18, LGALS9, and PTPN6) among immunologic genes, EMT genes, brown module genes, and module hub genes (Supplemental Fig. 4A). A heatmap illustrated the relationship between the expression of 23 IERGs, different subtypes, and clinical phenotypes, revealing significant positive associations of IS1, ES1, the socore of ESTIMATE, immune and stromal with these genes expression, while tumor purity showed a significant negative correlation (Supplemental Fig. 4B). Correlation analysis of the 23 IERGs demonstrated positive correlations among them (Supplemental Fig. 4C). KEGG pathway analysis identified the main enriched pathways (Supplemental Fig. 4D). The KEGG pathway classes were displayed using a lollipop plot (Supplemental Fig. 4E). A chord diagram illustrated the interaction between IERGs and KEGG enrichment (Supplemental Fig. 4F). Further functional analysis of the 23 IERGs using GO enrichment identified three principal categories: biological process (BP), cellular component (CC), and molecular function (MF) (Supplemental Figs. 4G–I). The BP category included defense response, cell activation, and immune response-related items (Supplemental Fig. 4G). The CC category featured secretory granules, secretory vesicles, and cell surface items (Supplemental Fig. 4H), while the MF category highlighted receptor regulator activity, signaling receptor binding, and cytokine receptor binding items (Supplemental Fig. 4I).

2.6 Prognostic analysis of 23 IERGs and a survival predictor model based on 16 IERGs

The analysis of univariate Cox regression was conducted for exploring the interactions between 23 IERGs expression and OS in glioma patients within the TCGA cohort, identifying 16 IERGs significantly associated with OS (Fig. 3A). To further refine this association, the optimal lambda value was calculated using LASSO-Cox regression with tenfold cross-validation based on the minimum partial likelihood deviance (Fig. 3B), highlighting the 16 IERGs linked to OS (Fig. 3C). Coefficient analysis revealed that CCR5, IRF8, EGF, FBP1, CXCR2, IL18, C5AR1, HAVCR2, PDGFB, CCL22, and CTSZ had positive coefficients, indicating their role as risk factors for gliomas, while IL6R, IL27, DAB2, NLRP3, and IGF1 had negative coefficients, suggesting their protective role (Fig. 3C). The 16-gene signature model was validated across the TCGA (Fig. 3D), REMBRANDT (Fig. 3G), and GRAVENDEEL (Fig. 3J) cohorts. In the TCGA (Fig. 3E), REMBRANDT (Fig. 3H), and GRAVENDEEL (Fig. 3K) cohorts, a significant correlation between higher risk scores and worse survival rates was found using Kaplan-Meier survival analysis. ROC curves assessed the accuracy of predictive model for 1-, 3-, and 5-year survival status, using data from the TCGA (Fig. 3F), REMBRANDT (Fig. 3I), and GRAVENDEEL (Fig. 3L) cohorts. AUC values for the TCGA cohort were 0.89 at 1 year (95 % CI, 0.87–0.92), 0.92 at 3 years (95 % CI, 0.89–0.95), and 0.84 at 5 years (95 % CI, 0.78–0.89), demonstrating high reliability of the 16 IERGs risk score as a prognostic indicator (Fig. 3F). Similarly, the REMBRANDT cohort showed favorable predictive AUC values of 0.71 at 1 year (95 % CI, 0.64–0.78), 0.8 at 3 years (95 % CI, 0.75–0.85), and 0.79 at 5 years (95 % CI, 0.74–0.85) (Fig. 3I). The GRAVENDEEL cohort presented AUC values of 0.82 at 1 year (95 % CI, 0.77–0.87), 0.86 at 3 years (95 % CI, 0.82–0.9), and 0.84 at 5 years (95 % CI, 0.78–0.89), further substantiating the 16 IERGs risk score's effectiveness as a prognostic indicator (Fig. 3L).Fig. 3 Prognostic analysis of 23 IERGs and a survival predictor model based on 16 IERGs. (A) Forest plot showing univariate Cox regression analysis results for the association between the expression of 23 IERGs and OS in patients with glioma from the TCGA cohort. (B) Least absolute shrinkage and selection operator (LASSO) analysis for the 23 IERGs, displaying the log (Lambda) values and identifying the optimal log (Lambda) value. (C) Coefficients of the 16-gene signature derived from the LASSO model. (D-F) Prognostic analysis of the 16-gene signature in the TCGA cohort: distribution of risk scores and gene expression levels (D), Kaplan-Meier survival curve comparing OS between high-risk and low-risk groups (E), and time-dependent ROC curve analysis for evaluating the 1, 3, and 5-year AUC of the risk score for OS prediction (F). (G-I) Prognostic validation of the 16-gene signature in the REMBRANDT cohort. (J-L) Prognostic validation of the 16-gene signature in the GRAVENDEEL cohort.

Fig. 3

2.7 Independent prognostic value and immune-related signature of 16-IERG risk group

A compound nomogram was established based on risk score and clinical characteristics to predict the 1-, 3-, and 5-year OS of gliomas patients, in order to create a more dependable predictive tool for clinical practice. (Fig. 4A). To evaluate the degree of agreement between anticipated and observed results at these time points, a calibration plot was made. (Fig. 4B). ssGSEA elucidated the relative immune infiltration levels within glioma samples from the TCGA cohort. The results revealed positive correlations between different immune score, IS1 and ES1 signatures, risk score, and immune cell infiltration (Fig. 4C). Greater immune cell infiltration was seen in high-risk locations for 16 IERGs. In the high-risk group, the abundance of anti-tumor and pro-tumor suppressive cells was positively correlated based on Pearson's correlation analysis (Fig. 4D), suggesting that anti-tumor inflammation may recruit or differentiate immunosuppressive cells. Comparison of tumor immune parameters between two distinct risk groups for 16 IERGs showed that the high-risk exhibited elevated TIDE scores (Fig. 4E), immune scores (Fig. 4G), stromal scores (Fig. 4H), and ESTIMATE scores (Fig. 4I), while demonstrating decreased stemness scores (Fig. 4F). This implies that glioma patients in the low-risk group might have more tumor differentiation, whereas those in the high-risk group might respond less well to ICB immunotherapy and have shorter post-treatment survival.Fig. 4 Independent prognostic value of the risk signature and ssGSEA for risk score in the TCGA cohort. (A) Nomograms for predicting the 1-year, 3-year, and 5-year survival probability of patient mortality based on the risk score and clinical variables, including age, gender and grade. (B) Plots illustrating the calibration of nomograms based on the risk score in terms of agreement between predicted and observed 1-year, 3-year, and 5-year outcomes. (C) ssGSEA analysis to identify the relative infiltration of immune cell populations in the TCGA cohort, with each cell type's infiltration normalized into a z-score. (D) Correlation between the infiltration of cell types involved in anti-tumor immunity and those performing pro-tumor, immune suppressive functions. (E-I) Comparisons of the TIDE score (E), stemness score (F), immune score (G), stromal score (H), and ESTIMATE score(I) in the high-risk and low-risk groups. The significance of the two groups was determined using the Wilcox test. P-values were showed as: ****, P < 0.0001.

Fig. 4

2.8 Single-cell TME characteristics

Single-cell TME characteristics of patients with glioma were investigated using the TISCH database, exploring the GSE131928 dataset comprising 13,552 cells based on the 10X Genomics platform. UMAP plots depicted cells colored by cluster (Fig. 5A), age (Fig. 5B), gender (Fig. 5C), and cell types (Fig. 5D–F). The initial analysis robustly grouped 13,552 glioma cells into 27 clusters (Fig. 5A), allowing intuitive observation of age, gender, and cell type distribution across these clusters (Fig. 5B–F). Pie and bar plots displayed the cell number of each type (Fig. 5G), as well as the cell-type proportion for each patient (Fig. 5H). Six immune checkpoints were examined for their expression levels across different cell types. CTLA4, HAVCR2, and PDCD1 were detected in CD8Tex cells; HAVCR2 existed in M1 cells; HAVCR2 and PDCD1LG2 were found in monocyte cells (Fig. 5I). The distribution of four hub IERGs (HAVCR2, IL18, LGALS9, and PTPN6) across different cell types indicated they were predominantly expressed in CD8Tex cells, M1 cells, and monocyte cells (Fig. 5J and K).Fig. 5 Single-cell tumor microenvironment characteristics. (A-H) UMAP plots representing cells colored by cluster (A), age (B), gender (C), and different cell types (D-F), and the pie chart (G) and bar plot (H) illustrating the cell number of each cell type and the cell-type proportion for each patient, respectively. (I) The expression levels of 6 immune checkpoints (CD274, CTLA4, HAVCR2, IDO1, PDCD1 and PDCD1LG2) in different cell types. (J-K) The expression distribution of 4 hub IERGs (HAVCR2, IL18, LGALS9 and PTPN6) in different cell types.

Fig. 5

2.9 Prognostic value of 4 hub IERGs in gliomas

In the TCGA cohort, the relationship between glioma patients OS and four hub IERGs expression was studied with univariate Cox regression analysis. The found revealed a relationship between OS and the expression of all four hub IERGs (Fig. 6A). We found that only IL18 expression was significantly related with OS based on a multivariate Cox regression analysis (Fig. 6B). These results indicate that IL18 can serve as a prognosticator in gliomas independently. Consequently, a predictive nomogram that included IL18 expression, age, gender, and grade was utilized to estimate survival probability at 1-, 3-, and 5-year intervals. (Fig. 6C). This nomogram disclosed strong agreement between observations and predictions using the calibration plot (Fig. 6D). Patients with glioma exhibiting elevated IL18 expression experienced significantly lower OS, demonstrating a negative association between IL18 expression and OS (Fig. 6E). Likewise, glioma patients with elevated IL18 expression showed a correlation with reduced progression-free survival (PFS), further emphasizing its negative association with PFS (Fig. 6F).Fig. 6 Prognostic value of 4 hub IERGs in gliomas. (A–B) Univariate (A) and multivariate (B) Cox regression of the 4 hub IERGs in gliomas. (C) Nomograms for predicting the 1-year, 3-year, and 5-year survival probability of patient mortality based on the results of multivariate Cox regression and clinical variables such as age, gender and grade. (D) Plots depicting the calibration of nomograms based on IL18 in terms of agreement between predicted and observed 1-year, 3-year, and 5-year outcomes. (E–F) Kaplan-Meier survival analysis of OS (E) and PFS (F) in the IL18 high-expression group and low-expression group.

Fig. 6

2.10 IL18 expression in gliomas

Using the TCGA and GTEx databases, we conducted an analysis of IL18 levels in gliomas and normal brain. In comparison to normal brain tissue, LGG and GBM exhibited significantly higher levels of IL18, with GBM exhibiting higher IL18 levels than LGG (Fig. 7A). To validate these results, we quantified the IHC results of IL18 protein expression in LGG and GBM samples from HPA database. Undetectable IL18 expression was observed in normal brain tissue and LGG samples, while GBM samples displayed medium IL18 staining (Fig. 7B). Further IHC staining revealed that Grade II gliomas exhibited lower IL18 protein levels compared to the markedly elevated levels found in Grade III and IV gliomas (Fig. 7C). Additionally, IL18 expression data from various glioma cell lines in the CCLE database were visualized (Fig. 7D). Validation of these findings through assessment of IL18 expression showed high IL18 level in U-251 MG cells and undetectable levels in HEB cells (Fig. 7E).Fig. 7 High expression of IL18 in GBM samples. (A) IL18 expression levels in LGG (n = 513), GBM (n = 153) and normal tissues (n = 2642) from the GTEx database. (B) IL18 expression in normal brain tissue (n = 3), LGG (n = 4) and GBM (n = 7) from HPA database. (C) IL18 expression in Grade II, Grade III and Grade IV glioma samples scale bar = 60 μm. (D) IL18 expression in different glioma cell lines from CCLE database. P-values were indicated as: *, P < 0.05; ****, P < 0.0001.

Fig. 7

3 Materials and methods

3.1 Data acquisition and procession

Sourced from The TCGA (http://cancergenome.nih.gov/abouttcga) and GlioVis (http://gliovis.bioinfo.cnio.es/) databases, we utilized three public glioma cohorts containing level 3 RNA-seq data (TPM-normalized expression matrix) along with corresponding clinical and follow-up information. Data from 663 TCGA-LGG/GBM samples were analyzed to identify prognostic gene expression signatures and construct prognostic risk signature models. The risk model was evaluated for its efficacy using the validation cohort from REMBRANDT dataset (n = 444) and GRAVENDEEL dataset (n = 276) provided by the GlioVis database. RNA-seq data were TPM-normalized to ensure cross-sample comparability. Clinical and follow-up information was meticulously curated to align with the RNA-seq data, including only samples with complete data for comprehensive analysis.

3.2 Unsupervised consensus clustering of gliomas

We retrieved 1959 immunologic genes and 1263 EMT genes from the ImmPort database (https://immport.niaid.nih.gov) and the dbEMT2.0 database (http://dbemt.bioinfo-minzhao.org/download.cgi). Novel immune and EMT classification subtypes for gliomas in the TCGA cohorts were delineated by unsupervised consensus clustering, using the immunologic and EMT genes with the R package “ConsensusClusterPlus v1.54.0” [18]. The maximum number of clusters was set to six, with 80 % of the total samples drawn 100 times (clusterAlg = “hc”, innerLinkage = ‘ward.D2′). Additionally, TCGA cohorts were analyzed for the distribution of tumor mutational burden (TMB) and microsatellite instability (MSI). TMB scores measure the number of mutations in tumor cells, with higher scores generally indicating a greater likelihood of response to immunotherapy due to a higher neoantigen load. MSI scores assess genetic hypermutability resulting from impaired DNA mismatch repair, with higher scores often correlating with better immunotherapy responses due to increased mutations and neoantigens. The relationships were shown using Spearman's correlation analysis, where a p-value of less than 0.05 was deemed statistically significant. Clinical information for the different ISs and ESs is provided in Supplemental Table 2.

3.3 Kaplan-Meier survival analysis

The OS of glioma patients was examined through Kaplan-Meier survival analysis, accompanied by a log-rank test, to evaluate differences based on subtypes and identified genes using the R packages “survival” and “survminer” (https://www.bioconductor.org/). P-values and 95 % confidence intervals (CIs) were calculated using log-rank tests and univariate Cox proportional hazards regression for Kaplan-Meier curves. R package “forestplot” was employed to generate the forest plot to illustrate the p-values and HR accompanied by 95 % CI for the identified genes (https://www.bioconductor.org/).

3.4 GESA and functional annotation

Based on p < 0.05 and FDR<0.25 as statistical basis, GSEA v3.0 (http://www.broadinstitute.org/gsea/) was used to conduct GSEA for HALLMARK terms in different subtypes based on expression profiles. The Molecular Signature Database v7.4 (http://www.broadinstitute.org/gsea/msigdb) was used to obtain the gene sets. With FDR <0.05 as the cutoff criterion, Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of identified genes were performed using “clusterProfiler” [19] package.

3.5 Prediction of immunotherapy responsiveness

Based on tumor transcriptome profiles in TCGA cohorts, the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was utilized to predict potential ICB responses [20]. The TIDE algorithm evaluates tumor immune evasion and its relationships with cytotoxic T lymphocyte infiltration levels in order to forecast patient survival and treatment response. It was developed to incorporate expression indicators of T cell dysfunction and T cell exclusion. Shorter survival following ICB treatment and a worsening curative effect are indicated by a higher TIDE score.

3.6 Evaluation of stemness index

The progressive loss of a differentiated phenotype and the acquisition of characteristics resembling stem cells and progenitors are hallmarks of cancer progression. In glioma samples, the mRNAsi score (stemness index) was assessed using the One-Class Logistic Regression (OCLR) algorithm [21,22]. Increased tumor dedifferentiation and biological processes present in cancer stem cells are linked to higher stemness indices. The 663 glioma samples were scored using the same Spearman correlation approach applied to the RNA expression data. By deducting the minimum value and dividing by a linear translation of the maximum value, the dryness index was mapped to the range [0,1].

3.7 Time landscape

The TCGA cohort's stromal and immune cell infiltration and tumor purity were assessed by using the Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm [23]. The ESTIMATE score, which is a composite score that reflects tumor purity; the stromal score, which indicates the abundance of stromal cells; the immune score, which indicates the abundance of immune cells infiltrating the tumor; and the tumor purity score, which directly indicates tumor purity, are the four scores that are computed by this algorithm. A higher ESTIMATE score signifies lower tumor purity. Greater infiltration of immune and stromal cells independently is indicated by higher immune and stromal scores. According to gene expression data, 22 immune cell types were quantified in tumor samples using the CIBERSORT algorithm [24] to provide a detailed picture of immune infiltration landscape in tumors. Furthermore, the association of 25 ICD modulators (e.g., ANXA1, IFNB1), 47 immune checkpoints (e.g., CTLA4, PD-L1), and 24 human leukocyte antigen family genes (e.g., HLA-A, HLA-E) from a previous study [25] with subtypes and immune microenvironment scores was examined.

3.8 Identification of IERGs

IERGs were identified by constructing co-expression networks, and gene modules closely related to immunoscores were recognized using “WGCNA” package [26]. An appropriate soft threshold power (β) was used to convert the gene similarity matrix created by Pearson's correlation analysis into a scale-free co-expression network. After converting this matrix into a topological overlap matrix (TOM), the dissimilarity matrix (1-TOM) between genes was computed. A higher TOM value denotes a higher level of connectedness. Genes with strong correlations were put in the same module. IERGs were found to be the intersections of strong immune-correlated gene modules, EMT genes, and immunologic genes. Additionally, the association of IERGs with clinical phenotype and immunoscores was examined.

3.9 Development and validation of prognostic risk signature model

Combined with a Cox regression model, the least absolute shrinkage and selection operator (LASSO) is used to reduce dimensionality and regularize biomarker screening [27]. In order to mitigate over-fitting and identify IERGs related to prognosis for feature selection, LASSO-Cox regression analysis with 1000 bootstrap samples and 10-fold cross-validation was carried out using the R package “Glmnet,” [28], ultimately establishing a prognostic risk signature model. The patients' median risk score was used to categorize them into low- and high-risk categories. These groups' survival differences were examined using a log-rank test and Kaplan-Meier survival analysis. Using the R package “survivalROC” (https://www.bioconductor.org/), the risk score's predictive accuracy was evaluated using time-dependent receiver operating characteristic (ROC) analysis. Age, sex, and tumor grade were incorporated in a prognostic nomogram that was created using the 'rms' R package (https://www.bioconductor.org/), and univariate and multivariate Cox regression analyses were performed to find and validate independent clinical prognostic markers. Significant differences in ICB response, stemness score, and immunoscore were assessed using the Wilcox test in the groups at high and low risk. The effectiveness of the prognostic risk signature model was further validated by calculating the risk scores of glioma patients from the GRAVENDEEL and REMBRANDT datasets.

3.10 Single-sample GESA

The 28 immune populations were measured using single-sample GESA (ssGSEA) using the R package “GSVA” [29]. Previously published studies provided gene markers for immune cells encompassing 782 genes for each type [30]. Enrichment scores in the ssGSEA using the R package “GSVA” represented the relative infiltration degree of each immune cell type [29]. Pearson's correlation analysis identified correlations between the infiltration of cell types executing anti-tumor immunity (Activated CD4 T cells, Activated CD8 T cells, Central memory CD4 T cells, Central memory CD8 T cells, Effector memory CD4 T cells, Effector memory CD8 T cells, Type 1 T helper cells, Type 17 T helper cells, Activated dendritic cells, CD56bright natural killer cells, Natural killer cells, Natural killer T cells) and those executing pro-tumor, immune-suppressive functions (Regulatory T cells, Type 2 T helper cells, CD56dim natural killer cells, Immature dendritic cells, Macrophages, MDSCs, Neutrophils, Plasmacytoid dendritic cells).

3.11 Single-cell TME analysis

The detailed cell-type annotation and exploration across various cancer types are provided by the single-cell RNA-seq database [Tumor Immune Single-cell Hub (TISCH)] [31]. The TISCH database investigates the single-cell TME characteristics of glioma patients, specifically the GSE131928 dataset comprising 13,552 cells from the 10X Genomics platform.

3.12 IL18 expression detection

Level 3 expression profiles based on RNA-sequencing were obtained from TCGA, together with matching clinical data for low-grade gliomas (LGG, n = 513), high-grade gliomas (GBM, n = 153), and paracancerous samples (n = 5). Expression data set was extracted from the GTEx database (https://portal.gtexportal.org/home/datasets) for normal brain tissue (n = 2642). HPA database provided immunohistochemistry (IHC) staining data for normal brain tissue (n = 3), LGG (n = 4), and GBM (n = 6), and we calculated the H-score.

3.13 Immunohistochemistry (IHC)

From the pathology department of Southern Medical University Shenzhen Hospital, thirteen paraffin-embedded glioma specimens, collected from 2019 to 2022, were acquired. The clinical research ethics committee of the hospital gave ethical permission. Following the manufacturer's instructions, 4 μm thick tissue slices were immunohistochemically stained using the PV-9000 Polymer Detection System (ZSGB-BIO, Beijing, China). Polyclonal anti-IL18 rabbit antibody (1:50; Affinity, Wuhan, China) that had been pre-diluted was incubated on the slices for a whole night at 4 °C. The Motic DSAssistant Lite system was utilized to take pictures after the slides had been counterstained with hematoxylin. Afterwards, we calculated H-scores for each sample.

3.14 CCLE database analysis, cell culture, and Western blotting verification

The CCLE dataset (https://portals.broadinstitute.org/ccle) provided the mRNA expression matrix of glioma cell lines. The following cells were cultivated in DMEM media (Thermo Scientific Hyclone, Beijing, China) supplemented with 10 % heat-inactivated fetal bovine serum (FBS, Sijiqing Biotech Corp., Hangzhou, China): HEB cells, U-87MG (RRID: CVCL_0022), U-251MG (RRID: CVCL_0021), and T98G (RRID: CVCL_0556). We prepared whole-cell extracts and electrotransferred the extracted lysates to polyvinylidene fluoride (PVDF) membranes after 10 % SDS-PAGE. IL-18 expression was probed via Western blotting using an IL-18 antibody (1:1000, Abclonal, Wuhan, China). Within the last three years, all cell lines were verified to be mycoplasma-free using STR profiling.

3.15 Statistical analysis

All analysis was conducted by R 4.3.3, and p < 0.05 was used to establish statistical significance.

4 Discussion

Gliomas are among of the most aggressive brain malignancies, with a dismal prognosis and limited therapeutic options [3]. Immunotherapy has shown effectiveness in various tumor treatments [32] and has emerged as a promising yet challenging option for gliomas [4,33,34], primarily due to the complexity of the TME. Characterizing TIME is crucial for designing successful antitumor immunotherapies [35,36]. This study explored the molecular subtypes of gliomas, classifying them into two ISs and two ESs, both displaying highly consistent clinical phenotypes and TIME characteristics. These findings were validated across different glioma cohorts, enhancing our understanding of TIME characteristics and providing implications for personalized immunotherapy in patients with glioma.

Zhang et al. [11] established a risk signature for glioma with six IRGs (CANX, HSPA1B, KLRC2, PSMC6, RFXAP, and TAP1). Their research demonstrated more CD4+ T cells existed in the low-risk group, while high-risk group showed elevated M1-type macrophages, M2-type macrophages, and CD8+ T cells in glioma cohorts [11]. Another study highlighted highly-expressed IRGs in GBM, including TREM1, GBP2, IFITM2, CIITA, and TYROBP, which are associated with worse prognosis and increased glioma grade [12]. IRGs that have been linked to poor prognosis include GRIA2 and MAP2, while BST2, B2M, and TRIM21 have been found to favorably correlate with OS by Huang et al. [13]. GRIA1, GRIA2, and MAP2 were strongly related with the infiltration degree of CD8+ T cells, neutrophils, and dendritic cells; B cells, macrophages, and dendritic cells were associated with BST2 and B2M; and both B cells and neutrophils were associated with TRIM21 [13]. Deng et al. [37] examined IRGs expression and patterns of immune infiltration, identifying CCL27, CXCL10, CCL20, RPR2, and SAA1 as high-risk factors, and SSTR1, SSTR6, HTR1A, MCHR2, RXFP3, and CPPM1 as low-risk factors [37]. These studies collectively underscore the pivotal role of IRGs in immune cell infiltration and glioma patient survival.

EMT-related gene induction, transcription regulation, and signaling significantly contribute to glioma malignancy, aiding prognosis assessment and understanding of EMT epigenetic modifications in glioma [14]. The EMT-related genes increase significantly in 3D cultured glioma models compared to those in two-dimensional culture, resulting in higher drug resistance [38]. Pilocytic astrocytomas and GBMs with vascular proliferation overexpress hub EMT factors like SLUG and TWIST [16]. Inhibition of EMT by silencing the Cullin-7 (CUL7) gene has been shown to reduce glioma malignant progression [39]. Similarly, the EMT process, migration and invasion in SHG44 and A172 glioma cells were markedly suppressed by COL3A1 knockdown [17].

Considering the synergistic involvement of IRGs and ERGs in forming the glioma TME, this study integrated ERGs and IRGs as IERGs to analyze their roles in predicting glioma behaviors. Sixteen IERGs were identified, with CCR5, IRF8, EGF, FBP1, CXCR2, IL18, C5AR1, HAVCR2, PDGFB, CCL22, and CTSZ acting as risk factors, while IL6R, IL27, DAB2, NLRP3, and IGF1 function as protective factors against gliomas.

Results from Cox regression analysis suggested that IL18 emerged as an independent clinical indicator for glioma patients. The active 18-kDa IL18 is derived from the proteolytic maturation of a 24-kDa precursor [40]. IL18 expression is implicated in various normal and tumor cells [41]. IL18 overexpression in bone marrow stromal cells (BMSCs) significantly inhibited glioma xenograft growth in rats, prolonging survival [42]. By improving T-cell infiltration and long-term anti-tumor immunity, mesenchymal stem cells (MSCs) overexpressing IL18 can prevent glioma growth and prolong survival in rats with gliomas., suggesting IL18 as a potential adoptive immunotherapy for malignant glioma [43]. Additionally, intraperitoneal rIL18 can delay subcutaneous glioma growth, though it has no apparent effects on gliomas in the brain, indicating the therapeutic potential of intratumoral rIL18 administration [44].

Contrarily, Yeh et al. reported that IL18 secreted by ECM-activated microglia generated from C6 gliomas can promote C6 glioma migration via the NO/cGMP pathway [45]. Consistent with Yeh's findings [45], our results demonstrated that patients with glioma exhibiting high IL18 levels have low OS, suggesting IL18 as an intrinsic glioma-promoting factor. Recently, disulfiram and ritonavir, drugs that may inhibit the maturation of pro–IL18, have been proposed for co-administration with temozolomide in GBM treatment trials [41,46].

5 Limitations

While there are several limitations that may impact the interpretation of the results., this work offers important insights into the TIME and suggests potential glioma biomarkers. The datasets used, primarily from TCGA and GlioVis databases, might not fully capture glioma heterogeneity across different populations. These findings need to be validated across a broader demographic spectrum and more cohorts need to be recruited. Additionally, cohort heterogeneity, including variations in glioma subtypes, disease severity, age, and underlying health conditions, might have influenced the outcomes. These factors can affect the expression levels of IERGs and the overall immune response. Stratifying participants based on these variables in future research could reveal more nuanced relationships between gene expression profiles and clinical characteristics. Technical variability in transcriptomic analyses and biomarker detection methods also poses challenges. Replicating the results using different technological platforms and conducting independent validation studies would help confirm the findings' reliability. Although significant associations between certain IERGs and glioma prognosis were identified, the functional roles and mechanistic contributions of these genes need further elucidation. Experimental studies focusing on the biological pathways involving these genes and their interactions with the TME and host immune system are essential to fully understand their roles in glioma progression and their therapeutic potential.

6 Conclusion

In summary, based on an integrated analysis of IRG and ERG sets in glioma patients, this study identified distinctive classification subtypes that can be applied to the development and implementation of effective immunotherapies in gliomas.Abbreviations list

Abbreviations	Meaning	
TIME	Tumor immune microenvironment	
ISs	Tmmune subtypes	
EMT	Epithelial-mesenchymal transition	
ESs	Epithelial-mesenchymal transition subtypes	
GBM	Glioblastoma	
OS	Overall survival	
TME	Tumor microenvironment	
IRGs	Immune-related genes	
ERGs	Epithelial-mesenchymal transition-related genes	
TIME	Tumor immune microenvironment	
CDX-LIPO	Liposomal honokiol and disulfiram/copper codelivery system	
NK	Natural killer	
TCGA	The Cancer Genome Atlas	
ISs	Immune subtypes	
ESs	EMT subtypes	
IERGs	Immune and EMT-related genes	
GESA	Gene set enrichment analysis	
ICB	Immune checkpoint blockade	
TMB	Tumor mutational burden	
MSI	Microsatellite instability	
ICPs	Immune checkpoint proteins	
ICD	Immunogenic cell death	
BP	Biological process	
CC	Cellular component	
MF	Molecular function	
PFS	Progression-free survival	
TMB	Tumor mutational burden	
MSI	Microsatellite instability	
HR	Hazard ratios	
CI	Confidence intervals	
GO	Gene Ontology	
KEGG	Kyoto Encyclopedia of Genes and Genomes	
TIDE	Tumor Immune Dysfunction and Exclusion	
OCLR	One-Class Logistic Regression	
mRNAsi score	Stemness index	
ESTIMATE	Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data	
HLA	Human leukocyte antigen	
TOM	Topological overlap matrix	
LASSO	Least absolute shrinkage and selection operator	
ROC	Receiver operating characteristic	
ssGSEA	Single-sample GESA	
MDSCs	Myeloid-derived suppressor cells	
TISCH	Tumor Immune Single-cell Hub	
IHC	Immunohistochemistry	
PVDF	Polyvinylidene fluoride	
LGG	Lower-grade glioma	
CUL7	Cullin-7	
BMSCs	Bone marrow stromal cells	
MSCs	Mesenchymal stem cells	

Conflict of interest disclosure

All authors confirm that they have no conflicts of interest.

Funding statement

The present study was supported by the 10.13039/501100001809 National Natural Science Foundation of China , grant number 81471279 and 81171138 , and Key discipline construction within the institution in 10.13039/501100013853 Shenzhen Hospital of Southern Medical University [grant number: zcxm-2023-xz-0003-07 ].

Ethics approval statement

This study was reviewed and approved by the Ethics Committee of Shenzhen Hospital, Southern Medical University with the approval case reference: [No. NYSZYYEC20230075], dated [November 18, 2023]. The study utilized archived paraffin tissue samples from the biobank in Shenzhen Hospital of Southern Medical University, which is a retrospective study without any intervention. The ethics committee exempted informed consent. All procedures were in accordance with the Declaration of Helsinki and the ethical principles of the medical community.

Patient consent statement

Not applicable.

Data availability statement

The datasets generated and/or analyzed during the current study are available in the Cancer Genome Atlas (TCGA) database [http://cancergenome.nih.gov/abouttcga] and GlioVis database [http://gliovis.bioinfo.cnio.es/]. Other data generated or analyzed during this study are included in this published article. If you have any questions, please contact the correspondence author (15wwlin@alumni.stu.edu.cn).

CRediT authorship contribution statement

Wen-wen Lin: Writing – original draft, Formal analysis, Data curation. Wei-jiang Zhao: Writing – review & editing, Writing – original draft. Guan-yong Ou: Writing – review & editing, Writing – original draft, Formal analysis, Data curation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following are the Supplementary data to this article:Multimedia component 1

Multimedia component 1

Multimedia component 2

Multimedia component 2

Multimedia component 3

Multimedia component 3

Supplemental Fig. 1 The study design.

Supplemental Fig. 1

Supplemental Fig. 2 TIME features of different subtypes in the TCGA cohort. (A-H) Comparisons of tumor purity and the immune, stromal, and ESTIMATE scores in the two ISs (A-D) and two ESs.(E-H). (I and K) Heatmaps of the 22 types of immune cells estimated by the CIBERSORT algorithm in the two ISs (I) and ESs (k), with different colors representing the distribution trend in different samples. (J and L) Comparisons of 22 types of immune cells in the two ISs (I) and ESs (J). The significance of the differences between the two groups was determined using the Wilcox test. P-values were showed as: ns, no significance; *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Supplemental Fig. 2

Supplemental Fig. 3 Association of ISs and ESs with ICPs, ICD modulators and HLA family in the TCGA cohort. Asterisks indicate the levels of significance for the differences between ISs or ESs. The significance of the differences between the two groups was assessed using the Wilcox test. P-values were displayed as: ns, no significance; *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Supplemental Fig. 3

Supplemental Fig. 4 Identification of IERGs and functional enrichment analysis. (A) Upset plot showing the intersection genes of immunologic genes, EMT genes, brown modules genes and hub genes. (B) Heatmap depicting the expression of 23 IERGs in different subtypes and clinical phenotypes. (C) Correlation heatmap of the 23 IERGs. (D) Bubble plot showcasing the top 10 enriched KEGG pathway for the 23 IERGs. (E) Lollipop plot displaying the KEGG pathway class, and (F) chord diagram showing the correlation between IERGs and KEGG pathways. (G-I) Bubble plots dispaying the top 10 enriched GO biological processes (G), cellular components (H) and molecular function (I) for 23 IERGs.

Supplemental Fig. 4

Acknowledgments

We appreciate Xin-yu Qiao changing the format of the manuscript.The authors express their gratitude to the anonymous reviewers whose insightful critiques and recommendations enabled us to improve our manuscript. We appreciate Bullet Edits Limited's assistance with the manuscript's language editing and proofreading.

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e36986.
==== Refs
References

1 Ostrom Q.T. Cioffi G. Gittleman H. Patil N. Waite K. Kruchko C. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2012-2016 Neuro Oncol. 21 Suppl 5 2019 v1 v100 10.1093/neuonc/noz150 31675094
2 Montoya M. Gallus M. Phyu S. Haegelin J. de Groot J. Okada H. A roadmap of CAR-T-cell therapy in glioblastoma: challenges and future perspectives Cells 13 9 2024 10.3390/cells13090726
3 Chen R. Smith-Cohn M. Cohen A.L. Colman H. Glioma subclassifications and their clinical significance Neurotherapeutics : the journal of the American Society for Experimental NeuroTherapeutics 14 2 2017 284 297 10.1007/s13311-017-0519-x 28281173
4 Xu S. Tang L. Li X. Fan F. Liu Z. Immunotherapy for glioma: current management and future application Cancer Lett. 476 2020 1 12 10.1016/j.canlet.2020.02.002 32044356
5 Hanahan D. Coussens L.M. Accessories to the crime: functions of cells recruited to the tumor microenvironment Cancer Cell 21 3 2012 309 322 10.1016/j.ccr.2012.02.022 22439926
6 Zhao W.J. Ou G.Y. Lin W.W. Integrative analysis of neuregulin family members-related tumor microenvironment for predicting the prognosis in gliomas Front. Immunol. 12 2021 682415 10.3389/fimmu.2021.682415
7 Lin W.W. Ou G.Y. Zhao W.J. Mutational profiling of low-grade gliomas identifies prognosis and immunotherapy-related biomarkers and tumour immune microenvironment characteristics J. Cell Mol. Med. 25 21 2021 10111 10125 10.1111/jcmm.16947 34597473
8 Quail D.F. Joyce J.A. Microenvironmental regulation of tumor progression and metastasis Nat. Med. 19 11 2013 1423 1437 10.1038/nm.3394 24202395
9 Wang X. Gao M. Ye J. Jiang Q. Yang Q. Zhang C. An immune gene-related five-lncRNA signature for to predict glioma prognosis Front. Genet. 11 2020 612037 10.3389/fgene.2020.612037
10 Zheng Z. Zhang J. Jiang J. He Y. Zhang W. Mo X. Remodeling tumor immune microenvironment (TIME) for glioma therapy using multi-targeting liposomal codelivery Journal for immunotherapy of cancer 8 2 2020 10.1136/jitc-2019-000207
11 Zhang M. Wang X. Chen X. Zhang Q. Hong J. Novel immune-related gene signature for risk stratification and prognosis of survival in lower-grade glioma Front. Genet. 11 2020 363 10.3389/fgene.2020.00363 32351547
12 Kong Y. Feng Z.C. Zhang Y.L. Liu X.F. Ma Y. Zhao Z.M. Identification of immune-related genes contributing to the development of glioblastoma using weighted gene Co-expression network analysis Front. Immunol. 11 2020 1281 10.3389/fimmu.2020.01281 32765489
13 Huang S. Song Z. Zhang T. He X. Huang K. Zhang Q. Identification of immune cell infiltration and immune-related genes in the tumor microenvironment of glioblastomas Front. Immunol. 11 2020 585034 10.3389/fimmu.2020.585034
14 Tao C. Huang K. Shi J. Hu Q. Li K. Zhu X. Genomics and prognosis analysis of epithelial-mesenchymal transition in glioma Front. Oncol. 10 2020 183 10.3389/fonc.2020.00183 32154177
15 Takashima Y. Kawaguchi A. Yamanaka R. Promising prognosis marker candidates on the status of epithelial-mesenchymal transition and glioma stem cells in glioblastoma Cells 8 11 2019 10.3390/cells8111312
16 Mäder L. Blank A.E. Capper D. Jansong J. Baumgarten P. Wirsik N.M. Pericytes/vessel-associated mural cells (VAMCs) are the major source of key epithelial-mesenchymal transition (EMT) factors SLUG and TWIST in human glioma Oncotarget 9 35 2018 24041 24053 10.18632/oncotarget.25275 29844871
17 Yin W. Zhu H. Tan J. Xin Z. Zhou Q. Cao Y. Identification of collagen genes related to immune infiltration and epithelial-mesenchymal transition in glioma Cancer Cell Int. 21 1 2021 276 10.1186/s12935-021-01982-0 34034744
18 Wilkerson M.D. Hayes D.N. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking Bioinformatics 26 12 2010 1572 1573 10.1093/bioinformatics/btq170 20427518
19 Yu G. Wang L.G. Han Y. He Q.Y. clusterProfiler: an R package for comparing biological themes among gene clusters OMICS A J. Integr. Biol. 16 5 2012 284 287 10.1089/omi.2011.0118
20 Jiang P. Gu S. Pan D. Fu J. Sahu A. Hu X. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response Nat. Med. 24 10 2018 1550 1558 10.1038/s41591-018-0136-1 30127393
21 Lian H. Han Y.P. Zhang Y.C. Zhao Y. Yan S. Li Q.F. Integrative analysis of gene expression and DNA methylation through one-class logistic regression machine learning identifies stemness features in medulloblastoma Mol. Oncol. 13 10 2019 2227 2245 10.1002/1878-0261.12557 31385424
22 Malta T.M. Sokolov A. Gentles A.J. Burzykowski T. Poisson L. Weinstein J.N. Machine learning identifies stemness features associated with oncogenic dedifferentiation Cell 173 2 2018 338 10.1016/j.cell.2018.03.034 54.e15 29625051
23 Yoshihara K. Shahmoradgoli M. Martínez E. Vegesna R. Kim H. Torres-Garcia W. Inferring tumour purity and stromal and immune cell admixture from expression data Nat. Commun. 4 2013 2612 10.1038/ncomms3612 24113773
24 Newman A.M. Liu C.L. Green M.R. Gentles A.J. Feng W. Xu Y. Robust enumeration of cell subsets from tissue expression profiles Nat. Methods 12 5 2015 453 457 10.1038/nmeth.3337 25822800
25 Huang X. Tang T. Zhang G. Liang T. Identification of tumor antigens and immune subtypes of cholangiocarcinoma for mRNA vaccine development Mol. Cancer 20 1 2021 50 10.1186/s12943-021-01342-6 33685460
26 Langfelder P. Horvath S. WGCNA: an R package for weighted correlation network analysis BMC Bioinf. 9 2008 559 10.1186/1471-2105-9-559
27 Tibshirani R. The lasso method for variable selection in the Cox model Stat. Med. 16 4 1997 385 395 10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3 9044528
28 Engebretsen S. Bohlin J. Statistical predictions with glmnet Clin. Epigenet. 11 1 2019 123 10.1186/s13148-019-0730-1
29 Hänzelmann S. Castelo R. Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data BMC Bioinf. 14 2013 7 10.1186/1471-2105-14-7
30 Charoentong P. Finotello F. Angelova M. Mayer C. Efremova M. Rieder D. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade Cell Rep. 18 1 2017 248 262 10.1016/j.celrep.2016.12.019 28052254
31 Sun D. Wang J. Han Y. Dong X. Ge J. Zheng R. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment Nucleic Acids Res. 49 2021 D1420 10.1093/nar/gkaa1020 d30 33179754
32 Riley R.S. June C.H. Langer R. Mitchell M.J. Delivery technologies for cancer immunotherapy Nat. Rev. Drug Discov. 18 3 2019 175 196 10.1038/s41573-018-0006-z 30622344
33 Medikonda R. Dunn G. Rahman M. Fecci P. Lim M. A review of glioblastoma immunotherapy Journal of neuro-oncology 151 1 2021 41 53 10.1007/s11060-020-03448-1 32253714
34 Wang S.S. Bandopadhayay P. Jenkins M.R. Towards immunotherapy for pediatric brain tumors Trends Immunol. 40 8 2019 748 761 10.1016/j.it.2019.05.009 31229353
35 Thorsson V. Gibbs D.L. Brown S.D. Wolf D. Bortone D.S. Ou Yang T.H. The immune landscape of cancer Immunity 48 4 2018 812 10.1016/j.immuni.2018.03.023 30.e14 29628290
36 Bindea G. Mlecnik B. Angell H.K. Galon J. The immune landscape of human tumors: implications for cancer immunotherapy OncoImmunology 3 1 2014 e27456 10.4161/onci.27456
37 Deng X. Lin D. Zhang X. Shen X. Yang Z. Yang L. Profiles of immune-related genes and immune cell infiltration in the tumor microenvironment of diffuse lower-grade gliomas J. Cell. Physiol. 235 10 2020 7321 7331 10.1002/jcp.29633 32162312
38 Wang X. Dai X. Zhang X. Ma C. Li X. Xu T. 3D bioprinted glioma cell-laden scaffolds enriching glioma stem cells via epithelial-mesenchymal transition J. Biomed. Mater. Res. 107 2 2019 383 391 10.1002/jbm.a.36549
39 Xu J. Zhang Z. Qian M. Wang S. Qiu W. Chen Z. Cullin-7 (CUL7) is overexpressed in glioma cells and promotes tumorigenesis via NF-κB activation J. Exp. Clin. Cancer Res. : CR 39 1 2020 59 10.1186/s13046-020-01553-7 32252802
40 Kim K.E. Song H. Hahm C. Yoon S.Y. Park S. Lee H.R. Expression of ADAM33 is a novel regulatory mechanism in IL-18-secreted process in gastric cancer J. Immunol. 182 6 2009 3548 3555 10.4049/jimmunol.0801695 19265133
41 Kast R.E. The role of interleukin-18 in glioblastoma pathology implies therapeutic potential of two old drugs-disulfiram and ritonavir Chin. J. Cancer 34 4 2015 161 165 10.1186/s40880-015-0010-1 25963312
42 Li C.H. Jiao B.H. [Effect of bone marrow stromal cells transfected with interleukin 18 on growth of intracranial glioma in rats] Ai zheng = Aizheng = Chinese journal of cancer 26 1 2007 38 43 17222365
43 Xu G. Jiang X.D. Xu Y. Zhang J. Huang F.H. Chen Z.Z. Adenoviral-mediated interleukin-18 expression in mesenchymal stem cells effectively suppresses the growth of glioma in rats Cell Biol. Int. 33 4 2009 466 474 10.1016/j.cellbi.2008.07.023 18725309
44 Kikuchi T. Akasaki Y. Joki T. Abe T. Kurimoto M. Ohno T. Antitumor activity of interleukin-18 on mouse glioma cells J. Immunother. 23 2 2000 184 189 10.1097/00002371-200003000-00002 10746544
45 Yeh W.L. Lu D.Y. Liou H.C. Fu W.M. A forward loop between glioma and microglia: glioma-derived extracellular matrix-activated microglia secrete IL-18 to enhance the migration of glioma cells J. Cell. Physiol. 227 2 2012 558 568 10.1002/jcp.22746 21442623
46 Sloand E.M. Maciejewski J. Kumar P. Kim S. Chaudhuri A. Young N. Protease inhibitors stimulate hematopoiesis and decrease apoptosis and ICE expression in CD34(+) cells Blood 96 8 2000 2735 2739 11023506
