2025 Volume 48 Issue 9 Pages 1319-1324
Glioblastoma (GBM) is a highly aggressive and lethal brain tumor with very poor prognosis despite recent progress in multimodal treatments. Within glioma tissue, various niche cells such as macrophages and neutrophils form a unique glioma immune microenvironment (GIME) by interacting with heterogenous cancer cells, and this has been implicated in disease progression and responsiveness to immunomodulatory therapies. This study explores novel potential prognostic markers associated with the GIME using integrated bioinformatics analyses, including single-cell RNA-sequencing (scRNA-seq), and spatial transcriptome (ST) datasets of clinical GBM specimens. We first identified 42 genes as being associated with poor prognosis in GBM from 5 different cohorts, GBM vs. nontumor tissue, grade IV vs. grade II gliomas, isocitrate dehydrogenase (IDH)-wild-type vs. IDH-mutant variants, mesenchymal vs. proneural and classical subtypes, and hazard ratio for overall survival. Among these, 32 genes were positively correlated with ESTIMATEScore, infiltration of various immune cell types, expression of known immune-related genes, and representative immune-associated biological signals. On scRNA-seq analysis, 7 genes were relatively concentrated in tumor-associated macrophages rather than in malignant cells. ST analysis revealed that Collagen beta(1-O)galactosyltransferase 1 (COLGALT1), Integrin subunit beta 2 (ITGB2), and Myosin light chain 12A (MYL12A) were distributed in the interface between the tumor and the peritumoral area, overlapping with the expression of representative immune-related genes. These findings support the potential of COLGALT1, ITGB2 and MYL12A as biomarkers for predicting the prognosis and immune responses of GBM, which can help in the development of potential immunotherapeutic strategies for GBM.
Gliomas are the most common type of primary central nervous system (CNS) cancer. WHO categorizes gliomas based on their level of malignancy: grades I and II are considered low-grade tumors, while grades III and IV are high-grade tumors.1) Glioblastoma (GBM), classified as a grade IV tumor, is a highly aggressive and lethal CNS cancer, with mesenchymal, classical, proneural, and neural subtypes.2) Isocitrate dehydrogenase (IDH) serves as a prognostic biomarker for patients with GBM.3) The IDH-mutant variant is typically associated with better overall survival rates, and responds better to temozolomide, a commonly used alkylating agent for both primary and recurrent high-grade gliomas. There have been recent advancements in multimodal standard treatment approaches with a combination of surgical resection, radiotherapy, and alkylating agent chemotherapy.4) However, patients with GBM still have very poor prognosis, with long-term survival being quite rare.
Immunotherapeutic strategies such as immune checkpoint inhibitors, adaptive T cell therapies, cancer vaccines, and oncolytic viral therapies have drastically improved outcomes for various types of cancers over the last decade, but these benefits are limited in gliomas, with most trials yielding negative results.5) Within the glioma tissue, heterogeneous niche cells interact with heterogenous cancer cells to form a unique tumor microenvironment (TME).6) Various native and infiltrating immune cells such as microglia, tumor-associated macrophages (TAMs), myeloid-derived suppressor cells, tumor-associated dendritic cells, tumor-associated neutrophils, and tumor-infiltrating lymphocytes all comprise a highly dynamic glioma immune microenvironment (GIME), which plays a crucial role in disease progression and responsiveness to immunomodulatory therapies.7)
Extensive research has sought to understand the pathophysiology of gliomas and determine reliable prognostic markers and therapeutic targets by utilizing bioinformatics technology. Bioinformatics analyses have paved the way for identifying potential prognostic markers and therapeutic targets in gliomas.8–10) Building on this knowledge, the present study explores the novel potential prognostic markers associated with GIME through integrated bioinformatics analyses, including bulk RNA-sequencing (bulk RNA-seq), single-cell RNA-seq (scRNA-seq), and spatial transcriptome (ST) datasets of clinical GBM specimens.
The bulk RNA-seq data for TCGA-GBM (The Cancer Genome Atlas Glioblastoma Multiforme) and TCGA-LGG (The Cancer Genome Atlas Low Grade Glioma) cohorts were downloaded using the GDCquery function in the TCGAbiolinks package (version 2.30.4). The CGGA (Chinese Glioma Genome Atlas) data (CGG_mRNAseq_693) were downloaded from the CGGA database (https://www.cgga.org.cn/index.jsp). The scRNA-seq data were obtained from GSE13504511) (Supplementary Fig. 1). The ST data were obtained from GSE25308012) (Supplementary Fig. 1).
Gene Expression ComparisonGene expression comparisons between cohorts were performed using TPM (Transcripts Per Million) and FPKM (Fragments Per Kilobase of exon per Million mapped fragments) values for TCGA-GBM and CGGA, respectively. Statistical significance was assessed using the Wilcoxon rank-sum test in the rstatix package (version 0.7.2). Multiple testing correction was performed using the Benjamini–Hochberg (BH) method.
Survival AnalysisPatients were divided into the high and low expression groups based on the median expression of each gene. Hazard ratios were calculated using the survival package (version 3.8.3), and the log-rank test was performed using the survdiff function from the same package.
TME EstimationFor the TCGA-GBM cohort, the ESTIMATE and ConsensusTME scores were calculated using the deconvolute_estimate and deconvolute_consensus_tme functions, respectively, in immunedeconv (version 2.1.0).13) Correlations between gene expression levels and these scores were computed using the corr.test function from the psych package (version 2.4.12).
Gene Set Enrichment Analysis (GSEA)GSEA was performed on the TCGA-GBM gene expression data after stratifying patients into high and low expression groups (based on median expression). The HALLMARK gene sets were downloaded using the msigdbr package (version 7.5.1). Differential expression analysis was conducted using DESeq2 on TCGA-GBM count data to calculate the log2 fold change (log2FC). Afterward, the genes were ranked based on log2FC in descending order. Normalized enrichment scores were computed using the GSEA function from the clusterProfiler package (version 4.10.1). Multiple testing correction was applied using the BH method.
scRNA-Seq AnalysisThe scRNA-seq data were analyzed using the Seurat package (version 5.2.1).14,15) After creating a Seurat object, cells with a mitochondrial gene content exceeding 20% were filtered out. The data were normalized using the SCTransform function, followed by principal component analysis. Integration was performed using the IntegrateLayers function with the PCAIntegration method. Clustering was conducted using the FindNeighbors, FindClusters (resolution = 0.5, algorithm = 1), and RunTSNE functions. Cell annotation was performed using the expression of marker genes for various cell types: malignant cells (SOX2, PTPRZ1, STMN2, MAP2, PTN), oligodendrocytes (TF, MAG, MOG, CLDN11, MBP, PLP1), endothelial cells (VWF, CDH5, ECSCR, SLCO2A1, MYCT1, TM4SF18), pericytes (RGS5), T cells (CD2, CD3D, CD3E, CD3G), and TAMs (CD14, AIF1, FCER1G, FCER3A, TYROBP, CSF1R). The distribution of gene expression was visualized using the FeaturePlot function.
ST AnalysisST data were analyzed using the Seurat package with patient GBM2-R1 data from GSE253080. After converting the data into Seurat objects, cells with mitochondrial gene content exceeding 30% or those with fewer than 1000 detected genes were removed. The data were normalized using the SCTransform function. Gene expression distributions were visualized using the SpatialFeaturePlot function with the data slot of the SCTAssay.
Statistical AnalysisStatistical analyses were performed using R software (ver. 4.3.0), with statistical significance mentioned accordingly in the results (* p < 0.05, ** p < 0.01, *** p < 0.001).
Differentially expressed genes (DEGs) in GBM were screened using TCGA and CGGA databases (Fig. 1A). There were 3,739 genes upregulated in GBM tissues (vs. nontumor brain tissues) in TCGA database, 3819 overlapping genes co-upregulated in grade IV tumors (vs. grade II) in TCGA and CGGA databases, 1665 overlapping genes co-upregulated in IDH-wild-type GBM (vs. IDH-mutant GBM) in TCGA and CGGA databases, and 2052 genes upregulated in mesenchymal GBM (vs. classical and proneural subtypes). Furthermore, 281 overlapping genes were identified among these 4 cohorts (GBM vs. nontumor tissue, grade IV vs. grade II gliomas, IDH-wild-type vs. IDH-mutant GBM, and mesenchymal vs. proneural and classical subtypes) (Fig. 1A). Univariate Cox regression analysis was conducted for these 281 genes to identify genes with potential prognostic value. Groups with high and low expression of these genes were compared, identifying genes with hazard ratio >1.0, resulting in 42 genes that were associated with poor prognosis in both the TCGA and CGGA cohorts (Figs. 1A, 1B).
(A) Venn diagram illustrating the overlap of DEGs identified across multiple comparisons, with 281 genes commonly identified in all comparisons: GBM tissues vs. nontumor brain tissues (TCGA, n = 296 vs. 5); grade IV vs. grade II gliomas (TCGA, n = 296 vs. 257; CGGA, n = 249 vs. 188), IDH-wild-type vs. IDH-mutant GBM (TCGA, n = 254 vs. 22; CGGA, n = 190 vs. 49); and mesenchymal subtype vs. classical and proneural subtypes (TCGA, n = 87 vs. 68 and 49). (B) Forest plot displaying the results of univariate Cox regression analysis for the 281 overlapping genes. Among these, 42 genes with a HR >1.0 were identified as potential prognostic markers that were associated with poor survival in both TCGA and CGGA cohorts.
Two algorithms (ESTIMATE and ConsensusTME) were used to determine the association between 42 candidate prognostic biomarkers and GIME in TCGA and CGGA databases. The ESTIMATE algorithm determines tumor purity (i.e., comprehensive ratio of the immune and stromal components in TME), while the ConsensusTME algorithm estimates the abundance of tumor-infiltrating immune cells.16,17) Furthermore, we also explored the correlation of the expression of candidate prognostic markers with the expression of immune-related genes and immune-associated biological processes and pathways via GSEA. Among 42 candidate prognostic markers, 32 genes were positively correlated with all four analyses related to the GIME: ESTIMATEScore (the sum of the infiltrating stromal and immune cell content) (Fig. 2A), infiltration of various immune cell types including macrophages and T cells (Fig. 2B), expression of known immune-related genes (Fig. 2C, Supplementary Fig. 2), and representative immune-associated biological processes and pathways (Fig. 2D).
(A, B) Correlation of the expression of the 42 candidate prognostic markers with the (A) ESTIMATE score and (B) ConsensusTME score in the TCGA cohort (n = 296). (C) Correlation between the expression of candidate prognostic markers and the expression of immune-related genes in the TCGA cohort (n = 296). (D) GSEA results for the HALLMARK immune-related gene sets in candidate prognostic markers high expression groups.
A combination of scRNA-seq and ST analyses of clinical GBM specimens were used to determine the distribution of 32 genes positively associated with GIME. Six clusters were identified through t-SNE (t-distributed Stochastic Neighbor Embedding) analysis based on genetic profiles (Fig. 3A). Canonical markers were used to annotate different cell types: TAM, malignant cells, oligodendrocytes, T cells, pericytes, and endothelial cells (Fig. 3B). Notably, representative immune-related genes such as CSF1R, HAVCR2, LGALS9, TGFB1, TGFBR1, and V-set Immunoregulatory Receptor (VSIR) were relatively distributed in TAMs compared to malignant cells (Fig. 3C). Out of 32 genes, Collagen beta(1-O)galactosyltransferase 1 (COLGALT1), Fc fragment of immunoglobulin G (IgG) receptor IIb (FCGR2B), Integrin subunit beta 2 (ITGB2), Leukocyte immunoglobulin-like receptor B3 (LILRB3), Myosin light chain 12A (MYL12A), Plasminogen activator urokinase receptor (PLAUR), and Solute carrier family 11 member 1 (SLC11A1) were mainly distributed in TAM compared to malignant cells (Fig. 3D). Kaplan–Meier analysis revealed that the expression of these 7 genes was correlated with poor prognosis in patients with GBM (Supplementary Figs. 3A, 3B). Afterward, the ST data of patients with GBM were used to complement the results of scRNA-seq (Fig. 3E). COLGALT1, ITGB2, MYL12A, and 6 immune-related genes were concentrated in certain areas adjacent to the glioma tumor, specifically at the interface between the tumor and the peritumoral area (Figs. 3F, 3G).
(A) t-SNE plot of scRNA-seq data from clinical GBM specimens, identifying six distinct cell clusters based on transcriptional profiles (7 patients; n = 25,612, TAMs; n = 8262, malignant cells; n = 14264, oligodendrocytes; n = 2058, T cells; n = 464, pericytes; n = 296, endothelial cells; n = 268). (B) Violin plot of canonical markers in each cell cluster. (C, D) Feature plot of representative (C) immune-related genes (CSF1R, HAVCR2, LGALS9, TGFB1, TGFBR1, and VSIR) and (D) candidate prognostic markers (COLGALT1, FCGR2B, ITGB2, LILRB3, MYL12A, PLAUR, and SLC11A1). (E) Schematic of the ST analysis in a patient with GBM (GBM2-R1; n = 3679). (F, G) Spatial feature plot of COLGALT1, ITGB2, MYL12A, and six immune-related genes.
COLGALT1 is one of two essential galactosyltransferases that initiate collagen glycosylation, which is a posttranslational modification.18) COLGALT1 has been reported as a potential biomarker for predicting the prognosis and immune responses of renal clear cell carcinoma.19) Meanwhile, MYL12A is a regulatory component in the myosin light chain family that controls myosin activity by phosphorylating the head region of myosin heavy chains. MYL12A has been associated with the development of breast cancer, because of its strongly positive expression in highly malignant invasive ductal carcinomas of the breast.20) However, the association of COLGALT1 or MYL12A with the pathogenesis of GBM remains elusive. Through integrated bioinformatics analyses, the present study has identified COLGALT1 and MYL12A as potential prognostic markers in GBM, due to their association with the GIME. Both were also closely related to immune cell infiltration, immune-associated signals, and the expression of immune-related genes. Furthermore, scRNA-seq analysis revealed that the expressions of COLGALT1 and MYL12A were concentrated in TAM rather than in malignant cells, while ST analysis revealed their distribution in the interface between the tumor and the peritumoral area, overlapping with the expression of representative immune-related genes.
ITGB2, the beta component of the β2 integrins, has been implicated in the pathogenesis of various cancers, such as liver cancer, colon cancer breast cancer, and leukemia.21) ITGB2 is also known as a predictor for the clinical prognosis and immunotherapy responsiveness of gliomas.22) Using ST and scRNA-seq analyses, the present study revealed that ITGB2 is distributed in interface between the tumor and the peritumoral area where TAM and T cells may be concentrated. This suggests that ITGB2 could also be a viable prognostic marker associated with GIME.
FCGR2B, a member of the fragment crystallizable receptor family, is recognized as a prognostic and immune microenvironmental marker in glioma patients. LILRB3, a member of the LILRB family, has been identified as an independent negative prognostic factor among various grades of gliomas. PLAUR, a component of the urokinase-type plasminogen activator system, functions as a regulator of the mesenchymal phenotype of GBM. SLC11A1, also known as natural resistance-associated macrophage protein-1, has been identified as a novel prognostic marker and as an indicator of immunotherapy response in gliomas. Although scRNA-seq analysis revealed that the expressions of FCGR2B, LILRB3, PLAUR, and SLC11A1 were concentrated in TAMs rather than in malignant cells, ST analysis did not show a clear distribution pattern for their gene expression.
This study has several limitations. We identified 281 upregulated genes among DEGs from four different cohorts (GBM vs. nontumor tissue, grade IV vs. grade II gliomas, IDH-wild-type vs. IDH-mutant GBM, and mesenchymal vs. proneural and classical subtypes). Although we focused primarily on upregulated genes due to their feasibility and translatability in developing inhibitors compared to activators, it should be noted that identifying downregulated genes could also provide valuable prognostic and mechanistic insights. We found that the distribution of the representative immune-related genes (CSF1R, HAVCR2, LGALS9, TGFB1, TGFBR1, and VSIR) was preferentially localized to TAMs rather than to malignant cells. Given that TAMs were the most abundant non-malignant cell type based on our analysis of scRNA-seq data from GSE135045, we focused our analysis on them in this study. However, further research is required to examine genes upregulated in immune cells other than TAMs.
Immunotherapeutic strategies have significantly improved the treatment of solid tumors in the past decade. These include immune checkpoint inhibitors (e.g., anti-programmed cell death-1 (PD-1)/PD-ligand 1 antibodies and anti-cytotoxic T lymphocyte-associated antigen-4 antibody) and chimeric antigen receptor-T cell therapy.23) Unfortunately, immunotherapy has limited benefits in GBM due to various reasons, such as the unique CNS immune profile, immunosuppressive TME, and heterogeneity of GBM cells.5) Through integrated bioinformatic analyses, this study aimed to identify pivotal prognostic markers of GBM and their association with the GIME, revealing the potential utility of COLGALT1 and MYL12A in immunotherapeutic strategies for GBM. Our findings carry substantial clinical implications for the treatment and prognosis of GBM, which can help in the development of immunotherapeutic strategies for the treatment of GBM and other malignant tumors associated with the TME.
This work was supported in part by the Japan Society for the Promotion of Science (20H03407 to E.H.).
The authors declare no conflict of interest.
This article contains supplementary materials.
