Bulk RNA-seq data, mutation data, and clinicopathological features for TCGA-GBM were downloaded from the UCSC Xena website (https://xena.ucsc.edu/). There are 168 glioma samples in the TCGA cohort. Gene expression profiling data for 374 GBM patients in the validation model were obtained from the China Glioma Genome Atlas (CGGA) data portal (http://www.cgga.org.cn/). All expression profile data were in TPM format. Batch correction and integration of the two sets of gene expression data were performed using the “limma” and “sva” (19) software packages. Single-cell RNA sequencing(scRNA-seq) data for GBM were downloaded from GSE84465 (20) and contained a total of 3589 cells within and near the tumor. Spatial transcriptome RNA sequencing(stRNA-seq) for primary GBM from the 10x Visium platform was downloaded from GSE194329 (21). 121 exosome-related genes were downloaded from the ExoBCD database (https://exobcd.liumwei.org) (22).

We conducted an analysis of single-cell RNA sequencing data using the R packages ‘Seurat’ and ‘SingleR’ (23). To ensure the inclusion of high-quality cellular data, we considered genes expressed in a minimum of three single cells. In addition, we excluded cells with gene counts fewer than 200 or exceeding 10,000, fewer than 1000 counts, and over 20% of mitochondrial and ribosomal genes. To address batch effects between cancer and paracancer samples, we employed the ‘harmony’ R package (24). The scRNA-seq data were normalized using the ‘Seurat’ R package with the ‘NormalizeData’ function. Subsequently, the normalized scRNA-seq data were transformed into Seurat objects, and the top 2000 highly variable genes were identified using the ‘FindVariableFeatures’ function. To reduce the dimensionality of the scRNA-seq data, we performed principal component analysis (PCA) using the ‘RunPCA’ function of the ‘Seurat’ R package. Significant principal components (PCs) were identified using JackStraw analysis, and appropriate PCs for cell clustering analysis were selected based on the proportion of variance. For the clustering of integrated data, we employed the ‘FindNeighbors’ and ‘FindClusters’ functions, and the resultant cells were visualized using UMAP or t-SNE methods. To identify genes specifically expressed in each cluster, we conducted Wilcoxon tests between pairs of cell clusters using the ‘FindAllMarkers’ and ‘FindMarkers’ functions of the ‘scran’ R package. The expression of specific genes was depicted using the ‘featureplot’ function (25). Cell type annotations were based on information from the original text and the tumor single-cell transcriptome database TISCH (http://tisch.comp-genomics.org/).

Cell developmental trajectories of malignant cells were analyzed by pseudotime using the “Monocle” R package (26). After transforming Seurat objects into cellular dataset objects, developmental difference genes were selected using unsupervised analysis. Utilizing the expression of 121 exosome-related genes, we used five commonly used algorithms to score gene sets from single-cell data (AddModuleScore, ssGSEA, AUCell, UCell, and singscore). The “AddModuleScore” algorithm, found within the “SingleR” R package, is used for scoring gene sets (27). Its essence lies in first calculating the mean of all genes in the gene set. Then, the expression matrix is divided into several segments based on the mean value, and a set of control genes is randomly sampled from each segment as background values. “ssGSEA” is a single-sample gene set enrichment analysis method used to evaluate the enrichment level of gene sets in a single sample or cell. It relies on the rank-based ordering of gene expression within a sample and calculates enrichment scores for each gene set. “AUCell” evaluates whether an input gene set is enriched among the top 5% expressed genes in a single sample, based on the ranking of gene expression in individual samples. The distribution of AUC scores across all cells enables the exploration of relative expression features. Given its rank-based scoring method, AUCell remains uninfluenced by gene expression units and normalization procedures. “UCell” (https://github.com/carmonalab/UCell) (28) is an unsupervised cell type identification method used to identify and classify the cell types of individual cells. UCell’s signature score is based on the Mann-Whitney U statistic, which is robust to dataset size and heterogeneity. Compared to other available methods, UCell requires less computational time and memory. “singscore” is a cell state assessment method used to quantify the activity level of specific functions or biological processes in a single sample or cell. It relies on gene sets from the gene expression profile and calculates cell state scores for samples or cells by considering gene weights and directions. These tools and methods have significant applications in single-cell transcriptomics and gene set enrichment analysis, helping researchers uncover changes in cell function, biological processes, and disease relevance, thus providing a deeper understanding of the complexity of biological systems (29).

The “CellCall” R package is a toolkit for inferring intercellular communication networks and internal regulatory signals by integrating intracellular and intercellular signals (30). The most notable feature lies in its ability to combine intercellular ligand-receptor communication with intracellular transcription factor expression, forming a ligand-receptor-transcription factor axis (L-R-TF axis). Simultaneously, it also encompasses pathway activity analysis, allowing for the analysis of receptor-cell pathway changes resulting from communication between two specific cell types.

The data is then analyzed using Seurat in R. The UMI counts are normalized, scaled, and the most variable features are determined by the “SCTransform” function. Dimensionality reduction is then performed using “RunPCA” for unsupervised clustering analysis. The “FindNeighbors” and “FindClusters” functions were executed with default parameters and the 30 most significant principal components. The “SpatialFeaturePlot” function was utilized for visualizing subgroups and genes.

The “scMetabolism” R package, developed by Fudan University, quantifies metabolic activity at the single-cell level. It operates based on a conventional single-cell matrix file and employs the VISION algorithm to score each cell, deriving a final activity score for each metabolic pathway (31).

For Python, the Scanpy and stlearn packages are employed. Scanpy is a Python-based package designed for analyzing single-cell data, encompassing pre-processing, visualization, clustering, proposed time series analysis, and differential expression analysis. The Institute of Molecular Biosciences at the University of Queensland has introduced an integrated analysis method, the stlearn package (https://github.com/BiomedicalMachineLearning/stLearn). This tool employs gene expression data, tissue morphology data, and spatial location information to initially identify cell types and subsequently reconstruct tissue cell types within tissues. It also infers evolutionary pathways and identifies tissue regions with high cell-to-cell interactions. stLearn integrates analyses to deduce interactions influenced by information on ligand pairs, gene expression, spatial location, and spatial cell type distribution.

Reverse Compositional Transcriptomics Deconvolution (RCTD) is a method utilized to deduce cellular composition in spatial transcriptomic data. This methodology involves analyzing the comprehensive gene expression profile of an entire tissue or sample to reverse-calculate the spatial distribution and relative abundance of each cell type. RCTD’s key advantage lies in its capacity to infer cellular composition from the overall gene expression profile, eliminating the necessity for the isolation and sequencing of individual cells. This attribute is particularly advantageous for spatial transcriptomics research, as it uncovers the relative distribution and interactions of distinct cell types within tissues.

The TCGA-GBM cohort was used as the training set, while the CGGA dataset was used as the validation set. First, 121 exosome-related genes were used as candidate predictors for univariate cox analysis to identify genes that were statistically associated with patient OS (p< 0.05). Next, we performed LASSO and multivariate regression analyses to further screen for genes and risk factors that were strongly associated with prognosis (32). The exosome-related index (ERI) was calculated for each GBM patient based on the coefficients determined by multivariate cox analysis. The risk score/ERI is calculated as follows:

Here, h0(t) represents the baseline hazard function, indicating the risk level when all predictor variables are at 0, serving as the increment over the risk at time t based on a zero value of predictor variables. b1, b2, …, bp are the regression coefficients of the Cox model, and x1, x2, …, xp are the corresponding predictor variables (gene expression values) at time t, the observed values. Based on the median value of ERI, the patients in the training and validation sets were divided into high and low-scoring groups. Survival curves were also plotted using the Kaplan-Meier method and log-rank tests were used to determine their statistical significance.

To calculate the probability of patient survival, we created a nomogram combining factors such as ERI, age, and IDH mutation status as prognostic factors. Using consistency index analysis and decision curve analysis (DCA), we further assessed the net benefit of the line graph and clinical characteristics alone. The ‘oncoplot’ functional waterfall plots from the R software ‘maftools’ package were used to explore detailed mutation characteristics. In addition, we investigated the correlation between ERI and tumor mutational load (TMB) and visualized this using the R software ‘ggplot2’ package (33).

Using the expression profile data, we employed the R package “estimate” to estimate the abundance of interstitial and immune cells, as well as tumor purity, in malignant tumor tissue (34). Subsequently, we assessed the extent of immune infiltration in GBM patients using the TIMER 2.0 database, which encompasses results from seven evaluation methods. These data were utilized to generate a heat map visualizing the relative fraction of immune cell infiltration within the tumor microenvironment (TME). To quantify the relative fraction of infiltrating immune cells and immune-related functions, we utilized the “ssGSEA” R package (35). Xu et al. created a website that provides us with a collection of genes related to cancer and immunology (36), as well as a collection of genes related to Mariathasan’s findings and a list of genes with favorable responses to anti-PD-L1 drugs (37). We utilized the GSVA approach to quantify both gene sets and evaluate their correlation with ERI. Furthermore, we explored the association between the model genes and the 51 immune genes, presenting the results in a circular heat map (38). To investigate variations in biological function across populations, we conducted Gene Set Enrichment Analysis (GSEA) using the MsigDB database “c2.cp.kegg.v6.2.” Additionally, we employed TIDE (http://tide.dfci.harvard.edu/), which stands for Tumor Immune Dysfunction and Rejection, as a computational framework to assess the potential for tumor immune escape based on the gene expression profile of tumor samples. The Cancer Immunome Atlas (TCIA) web tool provided comprehensive immunogenomic analysis results. The Immunophenotype Score (IPS), a quantitative tumor immunogenicity score ranging from 0 to 10, was used to denote the immunophenotype score (39). IPS can be used to predict response to immune checkpoint inhibitors.

U251MG, LN229, and SW1783 human glioma cells and human astrocytes (NHA) were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM, Gibco, C11995500BT, Canada) supplemented with 10% fetal bovine serum (FBS, Gibco, 10091148, Canada) and 1x penicillin/streptomycin (Gibco, 15140-122, Canada). All cultures were maintained in a CO2 incubator (TFS3111, USA) at 37°C with 5% CO2. BARD1 gene knockdown was achieved using small interfering RNA (siRNA). The specific BARD1 siRNA sequences can be found in Supplementary Table 1 . In brief, cells were seeded at 50% confluency in 6-well plates and transfected with negative control (NC) and siBARD1 using Lipofectamine 3000 (Invitrogen, USA).

Total RNA was extracted from cell lines and tissues using TRIzol (Sigma-Aldrich, T9424, America) according to the manufacturer’s instructions. cDNA was synthesized using the PrimeScriptTM RT Reagent Kit (Takara, RR047, Japan). Real-time polymerase chain reaction (RT-PCR) was performed using SYBR Green Master Mix (Q111-02, Vazyme) to quantify mRNA expression levels normalized to GAPDH mRNA levels. The 2−ΔΔCt method was used to calculate the expression levels. All primers were provided by Qingdao BioScience (Beijing, China), and the primer sequences can be found in Supplementary Table 1 .

First, cells (1000 cells per well) were seeded into a 96-well plate and incubated at 37°C for 4 hours with CCK-8 reagent (10 μL) (Dojindo, CK18, Japan). The absorbance was measured at a wavelength of 450 nanometers using an ELx800 plate reader (Thermo, Multiskan Spectrum, USA) to count the cells. Cell growth was represented as fold change from day 0 to day 4 and presented in a graph.

In a 6-well plate, transfected cells were electroplated and cultured in a cell culture incubator until 95% fusion. In each culture well, a straight line was scraped using a sterile 20 μL plastic pipette tip, and the unattached cells and debris were gently rinsed twice with PBS. Finally, the scratch wounds were photographed at 0 hours and 48 hours using Image J software, and the wound width was measured.

Cell invasion and migration studies were performed using a transwell assay. The upper chambers of a 24-well plate were filled with treated SW1783 cells (2×10^5 cells) and incubated for 48 hours. To evaluate the invasive and migratory abilities of the cells, the top surface of the plate was pre-coated with a matrix gel solution (BD Biosciences, USA) or left uncoated. The remaining cells at the bottom layer were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet (Solarbio, China) after removing the surface cells.

The statistical analyses were conducted using R version 4.1.3, 64-bit, along with its support packages. The Pycharm (3.9) integrated development environment for Python was also utilized. Prognostic values and comparisons of patient survival across different subgroups in each dataset were calculated using Kaplan-Meier survival analysis and the log-rank test. The non-parametric Wilcoxon rank sum test was employed to assess the relationship between two groups for continuous variables. Prognostic variables within the clinical characteristics of the different subgroups were identified through univariate and multivariate Cox regression analysis using the R package “survival.” Spearman correlation analysis was conducted to examine correlation coefficients. A significance level of P<0.05 was considered statistically significant for all statistical investigations.

A brief flowchart depicting this study is displayed in Figure 1 . Single-cell data from four samples were acquired based on the scRNA-seq data of GSE84465, aiming to explore the heterogeneity of glioblastoma and appraise disparities among cells situated within and adjacent to the tumor. Following stringent quality control filtration for the removal of low-quality cells, a total of 3533 cells were encompassed in subsequent analysis ( Supplementary Figure 1A ). After the elimination of batch effects and normalization of data, the top 2000 highly variable genes were chosen ( Supplementary Figure 1B ). For dimensionality reduction, the PCA technique was employed, with the top 20 principal components selected for further scrutiny based on P values ( Supplementary Figure 1C ). Data integration and batch effect elimination for the four single-cell samples were carried out using the Harmony algorithm, and t-SNE was employed to exhibit the distribution of cells before and after integration ( Supplementary Figures 1D, E ). Following the identification of 15 clusters, established marker genes were employed to annotate cell subpopulations. The distribution of different samples, tissue types, clusters, and cell subpopulations after annotation was demonstrated using t-SNE visualization ( Figures 2A–D ). Figure 2F delineates the proportional representation of distinct cell types within the four samples. The expression patterns of these genes across each cell type were showcased, predicated upon commonly used cellular marker genes ( Figure 2E ). The heatmap portrays the relative expression of marker genes within each cell subpopulation, with the top five marker genes calculated ( Figure 2G ). Grounded in the expression data of 121 exosome-related genes (ERGs), five prevalent algorithms (AddModuleScore, GSVA, AUCell, UCell, and singscore) were employed to assign scores to the gene sets within the single-cell data. The outcomes, depicted in Figures 2H, I , demonstrated relatively heightened exosome-related scores (ERS) in monocytes/macrophages and malignant cells. Moreover, a comparison between exosome-related scores in tumor and adjacent tissues unveiled intriguing observations: monocytes/macrophages, malignant cells, and oligodendroglial within tumors exhibited notably high ERS ( Figure 2J ). Following the integration of all malignant cells, a differential analysis of malignant cells within the tumor and adjacent tissues was conducted, revealing 33 ERGs with significant differential expression in both contexts ( Supplementary Table 2 ).

The evolution and differentiation of cells at the single-cell level can be inferred through cell trajectory analysis, accomplished by constructing intercellular trajectories to reshape cellular processes over time. The determination of cell trajectories and pseudotime distributions of malignant cells was carried out using the “monocle” R package, revealing the existence of seven cellular states during malignant cell development. For instance, the early state of cell development corresponds to cluster 5 ( Figure 3A ). The expression of 33 differentially expressed exosome-related genes at various developmental stages is depicted in the heatmap of Figure 3B , where BRAD1 is notably highly expressed in the early stages of malignant cell development. The “CellCall” tool was utilized to infer intercellular communication networks and internal regulatory signals by integrating intracellular and intercellular cues. By leveraging the median exosome-related scores across all malignant cells, we classified these cells into ERS high and ERS low groups. Bubble plots were employed to present the results of signaling pathway activity analysis, revealing that malignant cells with high scores exhibited intensified JAK-STAT signaling pathway activity in relation to monocytes/macrophages ( Figure 3C ). The JAK-STAT signaling pathway is implicated in pivotal biological processes like cell proliferation, differentiation, apoptosis, and immune regulation. The circle plot showcased in Figure 3D illustrates the intensity of direct ligand-receptor signaling across distinct cell types. Notably, a more pronounced CTF1-IL6ST ligand-receptor pair relationship was observed between malignant cells with high scores and monocytes/macrophages ( Figure 3F ). Subsequently, we shifted our focus to inferring the presence of ligand-receptor pairs and corresponding transcription factors (TFs) between high-scoring malignant cells and monocytes/macrophages ( Figure 3E ).

Spatial transcriptome sequencing data from a GBM patient was sourced from GSE194329, and following the exclusion of ribosomal and mitochondrial genes, the SCTransform method was applied to rectify sequencing depth and implement a series of normalizations. Through downscaled clustering, a total of 11 cellular subgroups were delineated within the spatial context ( Figures 4A, B ). The 33 exosome-related genes investigated in the prior study exhibited generally elevated expression across these 11 subgroups ( Figure 4C ). Guided by the original literature’s annotations, subpopulations 0 and 1 primarily occupied the GBM tumor core. Consequently, the “scMetabolism” R package was utilized for the assessment of metabolic activity within distinct cell subpopulations. Notably, subpopulations 0 and 1 exhibited close ties to Folate biosynthesis and tryptophan metabolic activity ( Figure 4D ). In this context, tryptophan, an essential amino acid, stands as a pivotal microenvironmental factor influencing the immunobiology of various tumor types. Its metabolism is integral to fostering immunosuppression, invasion promotion, and growth facilitation in malignant gliomas (40). As depicted in Figure 4E , heightened tryptophan metabolic activity was prominently enriched within the tumor’s core region. Leveraging Python’s Scanpy and stlearn packages, we conducted a trajectory analysis of cell subpopulations within the spatial domain. Employing pycharm, we normalized and clustered the spatial transcriptome data, yielding an additional set of 11 cell subpopulations ( Figure 4F ). Of interest, cluster 0, located within the tumor’s core, exhibited a differentiation trend toward cluster 5 situated in the periphery ( Figure 4G ). Ultimately, we employed the RCTD method to extrapolate the annotated cell types from single-cell data into the spatial dataset, inferring predominant cell types at each spatial location. Notably malignant cells with elevated exosome-related scores are predominantly located in the core of the tumor ( Figure 4H ).

To further utilize exosome-related features to aid clinical decision-making, 121 exosome-related genes (ERGs) were used to construct prognostic models in the bulk-seq of GBM. 168 GBM samples with survival information from the TCGA dataset were used as the training set for constructing prognostic risk models, and 388 GBM samples with survival information from the CGGA dataset were samples in the CGGA dataset used for external validation. First, the univariate cox analysis method was used to screen 17 ERGs affecting the overall survival (OS) of GBM patients (P<0.05, Figure 5A ). To avoid prognostic signature overfitting and narrowing down the genes predicting OS, a 6-gene signature was constructed by Lasso-Cox regression analysis ( Figures 5B, C ). Risk scores/ERI were calculated for each patient based on the expression profile of each gene and the corresponding regression coefficients ( Figure 5D ). By weighting the estimated cox regression coefficients, our model yielded an exosome-related index (ERI) for each patient in the TCGA cohort. Based on the scoring formula, patients were divided into a low-ERI group and a high-ERI group using the median ERI as the cut-off point. The risk graph shows the detailed survival outcomes for each patient in the TCGA cohort and the CGGA external validation cohort ( Figures 5E, F ). Survival curves for both the training and validation groups showed that patients in the high-ERI group had worse OS compared to those in the low-ERI group ( Figures 5G, H ). In addition, the ROC curves showed that ERI performed well in predicting OS for these individuals in the TCGA cohort ( Figure 5I ).

In the TCGA cohort, ERI could be an independent prognostic indicator for patients compared to other common clinical features (age, grade, IDH mutation status, etc.) according to univariate and multivariate Cox analysis ( Figures 6A, B ). In addition to that, the area under the curve (AUC) of ERI at 1 year was much higher than other clinicopathological features ( Figure 6E ). To help clinicians make better clinical decisions, based on the correlation between the above clinicopathological features and ERI, we created a nomogram for predicting 1-year, 1.5-year, and 2-year survival rates in GBM patients ( Figure 6C ). The calibration curves were also able to show that the nomogram was able to make accurate predictions ( Figure 6D ). 1-year DCA curves ( Figure 6F ) and C-index values ( Figure 6G ) both showed that our constructed nomogram and ERI had the highest net benefit and that the risk model constructed based on the six ERGs was more influential in clinical decision-making than the traditional model in clinical decision making. The results of the chi-square test showed that grouping was associated with the IDH mutation status of patients ( Figure 6H ). The proportion of IDH mutations was higher in patients in the low-ERI group ( Figure 6I ). Based on the analysis of the results, we are more confident that ERI and nomograms are reliable clinical predictive scoring systems.

Given the pivotal role of genetic mutations in tailoring cancer patient treatments, the somatic mutation profiles of 121 ERGs were examined. As depicted in Supplementary Figures 2A, B , HRNR emerged as the most frequently mutated gene, predominantly featuring missense mutations. Additionally, the distribution of the most frequently mutated genes in gliomas within ERI subgroups was scrutinized, revealing a higher prevalence of TP53 and EGFR mutations in high-ERI subgroups ( Supplementary Figure 2C ). Moreover, an exploration of the co-mutation patterns among model genes unveiled the co-mutation of HOXC6 and DCT (P<0.05, Supplementary Figure 2D ). However, no substantial discrepancy in TMB was observed between patients in the high-ERI and low-ERI groups ( Supplementary Figures 2E, F ). Subsequently, patients were categorized into four groups based on median TMB values and median ERI (H-TMB+high ERI, H-TMB+low ERI, L-TMB+high ERI, and L-TMB+low ERI). The outcomes indicated that patients with low ERI and high mutation burdens exhibited relatively improved prognoses ( Supplementary Figure 2G ).

The clinical outcome of patients and their response to therapy is influenced by the tumor microenvironment (TME). Among the factors within the TME, tumor-infiltrating immune cells play a significant role in impacting tumor progression and the effectiveness of antitumor therapies. Tumor-infiltrating immune cells (TIICs) constitute a crucial element of the TME, and their composition and distribution are intimately linked with tumorigenesis and progression (41). Consequently, an investigation was conducted into the immune landscape of high and low-ERI groups using algorithms from XCELL, TIMER, QUANTISEQ, MCPCOUNTER, CIBERSORT, CIBERSORT-ABS, and EPIC platforms ( Figure 7A ). To further delve into the relationship between the exosome-related index and immune cells and their functions, the enrichment scores of diverse immune cell subpopulations and immune functions were quantified through the “ssGSEA” method. The findings demonstrated that the high ERI group exhibited elevated scores for immune cell infiltration and immune pathway activation ( Figure 7C ). This encompassed certain immunosuppressive cells, including regulatory T cells (Tregs) and tolerogenic dendritic cells. Given the considerable influence of abnormal expression and function of immune checkpoint molecules on tumor immunotherapy, the expression of immune checkpoints (ICs) was evaluated across different ERI subgroups. Nearly all ICs showed heightened expression in the high ERI group ( Figure 7E ). Additionally, ESTIMATE was employed to compute the ratio of stromal and immune cells within various ERI subgroups, providing an estimate of tumor purity ( Figure 7F ). Heat maps were employed to visualize the expression of ICs, TME scores, and immune cell infiltration patterns in distinct subgroups ( Figure 7B ). Due to the contrasting prognostic outcomes and immune infiltration patterns in patients from the high-ERI and low-ERI groups, a Gene Set Enrichment Analysis (GSEA) was carried out to uncover potential disparities in biological functions between these two groups. For each group, the four most prominent signaling pathways were selected ( Figure 7D ). The high ERI group was notably associated with various cytokine-related pathways, whereas the low-ERI group exhibited enrichment in signaling pathways related to cell cycle and division. Based on the aforementioned findings, it is speculated that patients within the high ERI group may face a less favorable prognosis, although they exhibit more robust immune function. Patients with higher ERI might correspond to a tumor microenvironment in GBM characterized by immunosuppression, which could contribute to a reduced response rate to immunotherapy.

A recent study (42) has validated the remarkable advancements made in cancer treatment through immune checkpoint blockade (ICB) therapy. Gaining a deeper comprehension of the mechanisms underpinning cancer immunotherapy and translating this understanding into therapeutic strategies could potentially lead to prolonged survival for patients with limited treatment options. However, the efficacy of ICB therapy remains limited due to primary resistance, resistance development, and associated toxicities. To gain further insights into the role of risk scores in immunotherapy, the TIDE score was utilized to assess patients with potential abnormalities in immune function within tumors and regional lymph nodes. This approach aimed to judiciously identify candidates for immunotherapy. It is recognized that the impact of immunotherapy may fluctuate with tumor progression owing to variations in the degree of immune infiltration. Consequently, we probed whether prognostic models could predict responses to ICB in GBM patients. The utilization of the Immune Profile Score (IPS) can identify individuals who stand to benefit from immunotherapy. The violin plots illustrate the correlation between IPS values and risk groups, with higher IPS values indicating enhanced responses to PD-1 and CTLA-4 inhibitors. We anticipated a favorable immune response to CTLA-4 inhibitors in individuals from the low-ERI group ( Figure 8A ). Given the pivotal role of the immune microenvironment in mediating ICB responses, we further delved into the correlation between risk scores and ICB response characteristics. Notably, the ERI exhibited a significant negative correlation with these ICB response attributes ( Figure 8B ). In order to comprehensively explore the variations in immune responses across different subgroups, correlation analyses were conducted using six model genes and classical immune-related genes ( Figure 8C ). The TIDE score reflects tumor immune dysfunction and rejection, offering a computational framework to assess the likelihood of tumor immune evasion based on gene expression profiles from tumor samples. Elevated tumor TIDE prediction scores correlated with diminished ICB responses and inferior patient survival. Our findings revealed that patients in the high-ERI group exhibited heightened dysfunction scores, while those in the low-ERI group demonstrated elevated exclusion scores. Notably, there was no substantial discrepancy in TIDE scores between these two groups ( Figure 8D ).

Based on the intersection of differential and model genes in Supplementary Table 1 , a total of three ERGs are thought to influence GBM progression and treatment ( Figure 9A ). Interestingly BARD1, CTSB, and GSTP1 were all highly expressed in tumor samples from the TCGA dataset ( Figure 9B ). Based on whether these three genes were expressed in malignant cells, we divided the malignant cells in single-cell sequencing into expression positive and negative groups. Based on the marker genes of BARD1, CTSB, and GSTP1 expressing positive cells, the algorithm of GSVA was used to impute the content of these cells in the bulk sequencing data. We categorized them into high and low groups based on the optimal cutoff value, and the group with high expression of BARD1 and CTSB+ malignant cells had a poorer prognosis (P<0.05, Figure 9C and Supplementary Figure 3A ), whereas the difference in the overall survival curves between the two groups of GSTP1+ malignant cells was not statistically significant (P>0.05, Supplementary Figure 3D ). Patients with a lower percentage of BARD1+ malignant cells may respond better to immunotherapy (P<0.05, Figure 9D ). Meanwhile, the proportion of BARD1+ malignant cells was significantly positively correlated with cancer-associated fibroblasts (CAF) ( Figure 9E ). In contrast, patients with lower proportions of CTSB+ and GSTP1+ malignant cells may respond poorly to immunotherapy, and the comparison was not statistically significant (P>0.05, Supplementary Figures 3B, E ). In the spatial transcriptome, it is clear that BARD1 expression is overall low, while CTSB and GSTP1 are mainly expressed in the core region of the tumor ( Figure 9F , Supplementary Figures 3C, F ). Thus, malignant cells with positive BARD1 expression may be a risk factor for GBM. To understand the interactions of these cells in space, we again performed deconvolution of BARD1+ malignant cells, BARD1- malignant cells, and other cells according to the RCTD method. The recipient-ligand interactions were inferred by Python’s stlearn, and Figure 9G shows the top 50 reciprocal ligand-receptor pairs. Among them, there is strong ligand-receptor communication between BARD1+ malignant cells and vascular cells ( Figure 9H ). Figure 9I demonstrates the spatial enrichment score of COL4A2-CD93 ligand-receptor pairs.

The expression levels of BARD1 in HA cells and three GBM cell lines were initially compared through cell line experiments, revealing a notable upregulation of BARD1 in tumor cells ( Figure 10A ). Subsequently, the expression level of BARD1 was assessed 5 days post-transfection using qRT-PCR, to validate the effectiveness of siRNA-mediated BARD1 knockdown in the SW1783 cell line ( Figure 10B ). Following this, a CCK-8 cell assay demonstrated that the reduction of BARD1 resulting from knockdown significantly hindered the proliferative capacity of the SW1783 cell line ( Figures 10C, F ). Furthermore, the wound healing assay unveiled that the knockdown of the BARD1 gene notably impeded the migratory and invasive potential of SW1783 cells ( Figures 10D, G ). In alignment with the wound healing assay outcomes, GBM cells transfected with si-BARD1 exhibited diminished migratory and invasive capabilities in the transwell assay ( Figures 10E, H ). Collectively, these findings collectively indicate that BARD1 serves as a pro-carcinogenic factor in GBM.