Abstract This study presents a pioneering exploration into the role of aggrephagy-related genes (ARGs) in glioblastoma (GB), a kind of malignant tumor which is highly invasive and resistant to a series of therapy. Utilizing single-cell sequencing to dissect their influence on the tumor microenvironment (TME) and tumorigenesis. By applying non-negative matrix factorization for dimensionality reduction and clustering of single-cell data, distinct cellular subtypes within the TME influenced by ARGs were identified, uncovering their functions and interactions. The investigation extends to validating the prognostic significance of ARGs and their potential in predicting immunotherapy outcomes. Molecular docking analysis of key ARGs further highlights TUBA1C and UBB as promising therapeutic targets, offering novel insights into GB’s complex biology and suggesting a targeted approach for therapy, which is characterized by some crucial pathways in our analysis, including PI3k-akt and TGF-beta pathways. This comprehensive single-cell level examination not only advances our understanding of aggrephagy’s role in GB but also proposes new avenues for prognosis and treatment strategies, emphasizing the critical impact of ARGs on the TME and GB progression. Keywords: GB, Aggrephagy, scRNA-seq, Molecular docking, Drug target, Aggrephagy-related genes Introduction Glioma, a frequently encountered malignancy within brain tissue [[44]1]. GB, the most malignant among gliomas, is designated as grade 4 due to its distinctive histopathological features, encompassing necrosis and endothelial cell proliferation, as per the World Health Organization’s highest classification for brain tumors [[45]2]. Clinically, GB, being a high-grade glioma, for example, grade 4, remains a formidable challenge, characterized by its poor prognosis, and a 5-year survival rate of less than 10% [[46]3–[47]5]. Despite aggressive interventions including surgery, radiotherapy, and chemotherapy, the prognosis remains a persistent challenge [[48]6]. Hence, there exists a pressing need to delve deeply into GB’s pathogenesis and explore innovative cellular and molecular therapeutic approaches. At the cellular level, contemporary research increasingly acknowledges the pivotal role of the TME in fostering malignant tumor cell proliferation and advancement [[49]7–[50]9]. Within GB, the cells within the TME exhibit a remarkable degree of heterogeneity, coexisting symbiotically with tumor cells [[51]10]. Therefore, an exploration of GB’s evolution from a cellular vantage point within the TME can shed further light on the pathogenesis of GB and offer novel insights into GB immunotherapy and beyond. Differing from non-selective autophagy, aggrephagy is a selective process that discriminates and engulfs specific receptors, such as ubiquitin-tagged proteins [[52]10]. Aggrephagy especially targeted at the protein and promote cell life span compared with autophagy. In some cases, cells can trigger selective aggrephagy via NFκB, mTORC, and other pathways to eliminate aggregated proteins, playing a vital role in preserving cell survival and function [[53]11, [54]12]. Interestingly, aggrephagy may exert a positive regulatory influence on maintaining cellular homeostasis in tumor cells, potentially promoting tumor development. For instance, in hepatocellular carcinoma, aggrephagy is closely linked to various malignancies like HCC [[55]13]. Within the TME, aggrephagy is proposed as a crucial mechanism for upholding cellular homeostasis, thereby enhancing immunosuppression against tumor tissues [[56]14]. On one hand, aggrephagy aids in clearing misfolded proteins, thereby maintaining cellular homeostasis, particularly in tumor-associated macrophages (TAMs) that assist in immune evasion [[57]15]. On the other hand, during tumor angiogenesis, aggrephagy assists endothelial cells in maintaining equilibrium while executing essential functions like proliferation, differentiation, and protein synthesis, indirectly accelerating angiogenesis [[58]16, [59]17]. Several studies have already found that Aggrephagy plays a critical role in key genetic events for various cancers [[60]18]. Meanwhile, in GB, a highly malignant brain tumor, aggrephagy is closely intertwined with an essential pathway: the ubiquitin–proteasome system (UPS). This pathway can influence GB cell survival by regulating the metabolism and selective degradation of components related to cell death pathways [[61]19]. Nevertheless, in GB, the impact of aggrephagy on TME regulation and the interactions among various cell populations within the TME remain uncertain. Our work reveals the interaction by performing further investigations from cellular and molecular perspectives are warranted to elucidate the precise mechanisms through which ARGs influence GB and their potential implications for prognosis and immunotherapy. The aim of our study was to combine the analysis of scRNA-Seq in TME and bulk data in the TCGA-GB database to explore the mechanisms by which ARGs regulate the development of GB in TME at the single-cell level, as well as to investigate the impact of ARGs on the prognosis of patients with GB at the bulk level, and to further explore the GB patients’ immune response. The aim is to provide new perspectives and ideas for the prognosis of GB and to improve the prognosis of patient immunotherapy. Materials and methods Sources of data acquisition We utilized single-cell sequencing dataset from the [62]GSE162631 which is available in the GEO database ([63]www.ncbi.nlm.nih.gov/geo) for our single-cell analysis. Bulk-RNA data was obtained from The Cancer Genome Atlas (TCGA) databases ([64]https://portal.gdc.cancer.gov/), along with data from the GEO databases [65]GSE83300 and [66]GSE74187. These matrix data were integrated with clinical information. The collection of ARGs was acquired from the GSEA website ([67]https://www.gsea-msigdb.org/gsea/msigdb/cards/REACTOME_AGGREPHAGY) . Single-cell data reading and annotation We analyzed the tumor core and peritumor tissues from the dataset using the “Read10X” function within the R software’s “Seurat” package. Merging the two datasets with the “merge” function, we designated the single-cell clusters from the tumor core as tumoral clusters (TC) and the normal brain tissue samples from peritumor tissues as non-tumoral clusters (NTC). Following this, we conducted quality control by calculating the mitochondrial ribosomal gene ratio and erythrocyte ratio. Cells with gene numbers ranging from 200 to 4000 and mitochondrial ribosomal gene ratios below 15% passed quality control [[68]20]. Data normalization was achieved using the “NormalizeData” function in the “Seurat” package. Subsequently, we performed dimensionality reduction clustering by selecting the top 2000 highly variable genes with the “FindVariableFeatures” function. We set a principal component score from 1 to 10 for dimensionality reduction clustering, resulting in 20 clusters of cells. To annotate these clusters, we utilized marker genes from previous literature and data from the cellmarker website ([69]http://bio-bigdata.hrbmu.edu.cn/CellMarker/). Annotations were based on the expression of these genes in the 20 cell clusters. Cell communication CellChat analysis For cell populations of different species or subtypes, we used the CellChat package to reveal the strength of interactions and identify relationships between cell populations [[70]21]. We used the Seurat objects corresponding to the cell populations, created a CellChat object, and localized the ligand-receptors and signaling pathways playing important roles in the cell–cell interactions. Based on the ligand-receptor information in this database, we were able to localize the ligand-receptors and signaling pathways that play important roles in the interactions of each cell type with other cells. When there are fewer cell types generated by the classification, we also use the “netVisual_circle” function from CellChat to draw a circle diagram to demonstrate the counts and weights of cell interactions, while when there are too many subpopulations generated by the classification, we use the “netVisual_hierarchy1” function to draw a hierarchy plot to show the strength and number of interactions between different kinds of cells as receivers and targets, respectively. In order to further reveal the differences between the output and input signals of the interacting cells, we use the netAnalysis_signalingRole_heatmap function for intercellular signal comparison. Pseudotime analysis We employed the “Monocle2” package in R software for constructing pseudotime trajectories of individual cells. To explore the spatiotemporal expression changes in aggrephagy genes, we initiated the process by creating a mycds object as part of the Monocle2 package workflow using the “newCellDataSet” function within the Monocle package. We then employed the “dispersionTable” function to filter highly variable genes with average expression ≥ 0.1 and dispersion empirical > 1 * dispersion fit; thus, the number of highly variable genes is around 2000. Subsequently, we employed the DDRTree method from Monocle2 to dimensionalize the data, with a maximum of two components, and organized the data using the “orderCells” function. The outcomes of this analysis were visually represented using the “plot_pseudotime_heatmap” function within the Monocle package. In order to demonstrate their differences in pseudotime analysis more clearly and mapped the result with previous classification, we delve into the pseudotime trajectory of the cell population, projecting it onto a UMAP plot. This was achieved by merging the UMAP_1 and UMAP_2 columns from the downscaled Seurat object with the pseudotime data from the mycds object using the “FetchData” function in the Seurat package, facilitating the creation of a combined pseudotime trajectory and UMAP map. Single-cell differential gene analysis To discern differential gene expression between two cell populations, we harnessed the “FindMarkers” function within the Seurat package, which can find markers (differentially expressed genes) for identity classes. The method parameter was configured as “wilcox”. The p value was set below 0.05 as significant. Subsequently, we selected the top 10 differential genes to construct a heatmap and further visualized specific differentially expressed genes, including those associated with ARGs, using the “FeaturePlot” function, it can demonstrate 10 differentially expressed genes for either cluster. NMF for ARGs in TME cells in GB Non-negative Matrix Factorization (NMF) is a popular data dimension reduction method in recent years. The traditional NMF method has high sensitivity to data noise [[71]22]. Initially, we isolated the cell populations necessitating further subtype characterization and extracted their corresponding expression matrices. Subsequently, employing the NMF package in R software, we subjected these matrices to factorization, employing the NMF factorization results for dimensionality reduction clustering. And we visualized the partitioned NMF cell populations using the UMAP dimensionality reduction technique. Next, we utilized the “FindAllMarkers” function in the Seurat package, employing a log2FC threshold of 0.5 and an adjusted p value <0.05 to pinpoint the characteristic genes for each NMF cell cluster. Cell clusters were categorized based on these feature genes. If the first feature gene of a cluster belonged to ARGs with a log2FC greater than 1, it was designated as a metagene+ subtype [[72]23]. If the metagene did not pertain to ARGs, the first feature gene cluster was identified as a Non-aggregation subtype for ARGs. In cases where the metagene was ARG-related but the log2FC was less than 1, the cell was categorized as an Unclear subtype. Assessment of cellular metabolic activity We isolated a specific cell group from the Seurat object to assess cellular metabolism. Initially, we gauged the metabolic pathway strength across various cell subtypes using the “sc.metabolism.Seurat” function within the R software’s “scMetabolism” package [[73]24]. Pathways were annotated using KEGG entries. We opted to visualize the top 30 most active metabolic pathways, which were subsequently presented using the “DotPlot.metabolism” function [[74]25]. Single-cell GO enrichment analysis During the functional enrichment analysis of macrophage NMF subtypes, we harnessed the “ClusterGVis” package in R software to conduct GO enrichment analysis for macrophages across different NMF subtypes [[75]26]. We extracted the annotated macrophage population of NMF cell types and selected signature genes with log2fc > 0.25 and min.pct > 0.25, selecting the top ten marker genes with the highest expression [[76]27]. Data preparation was carried out using the “prepareDataFromscRNA” function, which calculates the average value for cells with the same gene. Subsequently, the “enrichCluster” function was applied to analyze the GO enrichment of cell clusters. The results of this enrichment were visually represented using the “visCluster” function, with parameters specifying plot.type = ”both,” add.bar = “T,” and the inclusion of “annoTerm.data” details. The visualization incorporated a combination of line plots, heatmap plots, annotation entries, and bar plots. For the functional enrichment of endothelial cells, we utilized the “irGESA” package in R software to perform GSEA on NMF cell clusters, opting for the UCell method from the available enrichment methods [[77]28]. Densitograms, half-violin plots, and other visual representations highlight gene expression variations in relation to specific pathways within NMF cell types. Single-cell regulatory network inference and clustering (SCENIC) analysis To construct single-cell gene regulatory networks at the single-cell level, we employed the SCENIC package in R software to identify transcription factors influencing cellular diversity among distinct subtypes [[78]29]. We imported two specific gene sequences from the RcisTarget database: hg19-500bp-upstream-7species.mc9nr and hg19-tss-centered-10kb-7species.mc9nr. Leveraging this sequence information, we independently delineated gene regulatory networks comprising transcription factors (TFs) and transcription start sites (TSS) for each cell subtype. The networks were visually represented through heatmaps [[79]30, [80]31]. Cell communication pathway identification and annotation Once we identified all NMF cell clusters comprising immune cells, we proceeded with cellular pathway enrichment analysis for each cell subtype. Our attention was particularly drawn to subtypes showing significant activity in intercellular communication in the enrichment results, warranting further investigation. To identify and analyze the key pathways contributing to intercellular communication through which these cell subtypes interact with others, we leveraged the “CommPath” package in R software. The initial step involved identifying ligand-receptor pairs crucial for cellular communication, and annotating pathways associated with the identified ligands using the KEGG database. Subsequently, we utilized the gsva method to compute the activation score for each pathway. Based on these scores, we pinpointed the pathways significantly up-regulated in each cell cluster, selecting the top ten for heatmap visualization. Upon reviewing the heatmap outcomes, we singled out the most active cell subtype in communication, identified as the central hub for communication. This central cell became the focal point for exploring essential cellular communication pathways it participated in. Our overarching objective was to build a comprehensive intercellular communication map, following this format: upstream cell-receptor-center cell-ligand-downstream cell. We also mapped the ligand-receptor interactions contributing to these pathways. ARGs regulated cell clusters score in Bulk RNA-sequence Having pinpointed cell subpopulations associated with ARGs within single-cell data, the prognostic and immunotherapeutic implications of these cell subpopulations warrant validation. In pursuit of this, we integrated TCGA-GB datasets into our analysis. We employed the “FindAllmarker” function in Seurat to identify marker genes for each NMF cell cluster and applied the GSVA (Gene Set Variation Analysis) to score the bulk data. Furthermore, we incorporated the clinical data corresponding to TCGA GB, stratified the samples into normal and tumor groups, and represented the scores for both in ARGs-regulated clusters. We depicted the abundance of each cell type in the tumor and normal groups using box plots. Survival analysis and K–M curves We excluded normal samples from TCGA and integrated survival time and status from clinical data with the TCGA GB matrix. Utilizing the marker genes of NMF cell clusters as the prognostic gene set, we employed the “survival” and “survminer” packages in R software for survival analysis. For each NMF cell cluster, we determined the group with high scores of the respective marker gene by finding the best cutoff when the gene scores surpassed the predefined threshold. Subsequently, we compared the impact of high scores versus low scores on the overall survival of GB patients and presented the results through K–M survival curves. In parallel, we employed GEO bulk data to replicate the same analysis, providing validation via K–M curves. Immune checkpoint blockade (ICB) Tumor Immune Dysfunction and Exclusion (TIDE) serves as a predictive tool for individualized immunotherapy response assessment in public cancer cohort [[81]32]. TCGA and GEO matrix data were independently submitted to the TIDE platform. The samples were classified as either responders or non-responders based on TIDE predictions and then visualized through bar plots. This approach effectively visualized the impact of NMF cell clusters on immunotherapy, using responder and non-responder scores to indirectly assess each cluster’s influence. The findings were presented via box plots. Molecular docking analysis Molecular docking analysis is a tool to demonstrate the interaction of 2 molecules. For example, a protein and a pharmaceutical small molecule. To assess the binding energies and interaction profiles between the drug candidates and our designated protein targets, we employed Vina 1.2.2, a computerized protein–ligand docking software. The molecular structure of Quercetin was sourced from the PubChem Compound Database [[82]33]. For proteins TUBA1C (PDB number, 8IXD) and UBB (PDB number, 2KX0), we acquired their 3D structures from the Protein Data Bank (PDB). Initial preparations involved converting all protein and molecule files to the PDBQT format, eliminating water molecules, and introducing polar hydrogen atoms. Grid boxes were meticulously aligned to encompass the structural domains of each protein while allowing for unhindered molecular mobility. The grid box dimensions were set to 30 Å × 30 Å × 30 Å cubic space with a grid point separation of 0.05 nm by AutoDock Tools. Subsequently, we conducted molecular docking investigations using Autodock Vina 1.2.2 ([83]http://autodock.scripps.edu/) for model visualization. Statistic analysis We conducted our routine analysis for this study using R 4.3.1 software. Our standard tests encompassed Student’s t test, Wilcoxon rank sum test, K–W test, and chi-square test. Pearson correlation was employed to ascertain the relationship between the expression of cells from various ARGs subtypes. A two-tailed p value <0.05 signified statistical significance. To delve into the linkage between critical TAM functions in GB and gene expression, we interrogated well-documented genes expressed by TAMs, along with a list of genes associated with TAM functions. Subsequently, we employed the “ComplexHeatmap” or “pheatmap” package to visualize our heatmaps, and we utilized violin and box plots for comparative analysis of data distribution between the two sets. Result The landscape of Aggrephagy in TME cells in GB The research methodology, illustrated in Fig. [84]1, involved the analysis of two samples from a GB patient to investigate the pivotal role of Aggrephagy in GB. One sample originated from the patient’s tumor core, labeled as TC, while the other came from peritumoral tissue, labeled as NTC. Following quality-controlled scRNA-seq, 26,669 cells were identified, and the primary constituents were selected for PCA down-clustering (Fig. [85]2A, B), resulting in 20 clusters (Fig. [86]2C). Comprehensive annotation revealed five distinct cell types, including macrophages, microglia cells, endothelial cells, monocytes, and T cells (Fig. [87]2D). We then examined intercellular communication, which was concentrated among macrophages, monocytes, and endothelial cells, while communication between T cells and microglia was relatively weak (Fig. [88]2E). Subsequently, we analyzed differentially expressed genes (DEGs) in GB samples, highlighting ARGs such as HSP90AA1, DYNLL1, UBC, and VIM, which exhibited distinct distributions in various subpopulations and high expression levels (Fig. [89]2F, G). Furthermore, we assessed the distribution of cell proportions in tumor and non-tumor samples. The results indicated that macrophages predominated in tumor samples, while microglia were more prominent in non-tumor samples. Additionally, the number and proportion of endothelial cells in tumor samples exceeded those in non-tumor samples (Fig. [90]2H). Fig. 1. [91]Fig. 1 [92]Open in a new tab Flow chart of the study Fig. 2. [93]Fig. 2 [94]Open in a new tab Cellular annotation based on ARG expression is illustrated. A The PCA dot plot distinguishes TC and NTC through various colors. B Principal component analysis was conducted, selecting 1:10 as the principal component score in the curve. C Dimensionality reduction clustering of single-cell data was achieved using the UMAP method for visualization. D The single-cell data were meticulously annotated, resulting in the classification of five distinct cell groups. E The communication status among these 5 annotated cell clusters is depicted. F The heatmap portrays ARG expression across these 5 cell clusters. H The percentage distribution of these 5 cell groups in NTC and the TC is presented. G Expression variations of select ARGs across the 5 cell populations are visually depicted using violin plots The alteration of phagocytes from normal brain tissue to tumoral brain tissue Indeed, microglia, the brain’s resident macrophages, are recognized for their trophic and phagocytic functions. Remarkably, in our patient, non-tumor tissues were predominantly inhabited by microglial cells, while macrophages dominated in tumor tissues. (Fig. [95]3A). Consistent with existing literature, one study reported macrophages constituting 85% of phagocytes in the TME, with microglia accounting for the remaining 15% [[96]34]. To gain insights into the dynamic changes within these phagocytes, we conducted pseudotime analysis. This revealed a progressive up-regulation of genes like UBA52, TUBB4A, TUBA4A, and DYNC1H1 during Aggrephagy-mediated cell development. Conversely, genes including TUBA1B, TUBB1, and ARL13B exhibited significant down-regulation over time (Fig. [97]3B). UMAP mapping clearly demonstrated that macrophage infiltration followed microglia, indicating a shift from microglia to macrophages in the TME (Fig. [98]3C). Further analysis of differentially expressed genes (DEGs) in the two types of phagocytes from the same patient uncovered notable distinctions. In macrophages, genes such as HSPA6, HSPA1A, SLC1A3, TSC22D3, and CACN1A exhibited significant down-regulation, while S100A family proteins were notably up-regulated (Fig. [99]3D, E). Additionally, several ARGs, including UBC, HSP90AA1, VIM, and TUBA1C, demonstrated significant differential expression (Fig. [100]3E). Fig. 3. [101]Fig. 3 [102]Open in a new tab Analysis of phagocytes in the tumor core and peritumor tissue revealed the following insights. A The distribution of infiltrated phagocytes differed between peritumor tissue and the tumor core. B Pseudotime analysis of macrophages and microglia was conducted to generate a heatmap of ARGs expression. C Pseudotime trajectories elucidated the chronological order of these two cell populations in the context of tumorigenesis, highlighting that macrophages appeared subsequent to microglia. D Differential gene analysis of macrophages and microglia involved the examination of ten up-regulated and down-regulated pathways in macrophages, verifying their immunosuppressive character. E Dot plots portrayed the distribution of genes, providing visual representation of the differing presence of certain ARGs and genes related to tumor immunity in these two cell populations Effects of macrophages regulated by ARGs in GB A total of 14,496 macrophages were identified in the scRNA-seq dataset. Following NMF dimensionality reduction clustering, eight distinct clusters emerged. These included the HSP90AA1-positive cluster, DYNLL1 cluster, TUBA1C cluster, PARK7 cluster, TUBB4B cluster, UBB cluster, an Unclear cluster, and the Non-Aggre cluster (Fig. [103]4A). To elucidate their metabolic characteristics, we delved into the metabolic pathway profiles. Metabolic pathway enrichment analysis revealed that the HSP90AA1-positive clusters, TUBA1C-positive clusters, PARK7 clusters, and unclear clusters exhibited more pronounced enrichment in metabolic pathways compared to the Non-Aggre clusters. These pathways encompassed fatty acid biosynthesis, steroid biosynthesis, sphingolipid metabolism, and glycerophospholipid metabolism (Fig. [104]4B). Further, we conducted functional enrichment analysis for these macrophage clusters, adding GO annotations to each. The results highlighted that HSP90AA1+mac, TUBA1C+mac, TUBB4B+mac, and UBB+mac exhibited elevated expression of functionally related genes, distinctly differing from Non-Aggre. Specifically, HSP90AA1+mac was primarily associated with cellular response to heat and unfolded protein, while TUBA1C+mac was mainly linked to the chemokine-mediated signaling pathway and neutrophil migration (Fig. [105]4C). Moreover, we conducted M1 and M2 typing of macrophages, revealing that Aggrephagy-positive clusters exhibited higher scores for M1 macrophages, the pro-inflammatory subtype (Fig. [106]4D). This suggests that aggrephagy may play a pivotal role in influencing the development and progression of GB by activating pro-inflammatory macrophage activities. Fig. 4. [107]Fig. 4 [108]Open in a new tab The exploration of macrophage metabolism and its regulation by ARGs led to the following findings. A Utilizing ARGs, macrophages were subjected to NMF-based subtyping, resulting in the categorization of ARGs subtypes. Visualization was achieved through the UMAP method. B Metabolic analysis of these ARGs subtypes in macrophages revealed variations in metabolic intensity, with point size indicating p values and point color representing metabolic intensity. C Functional enrichment analysis further illustrated differences in the expression of various gene modules across different subpopulations. Heatmaps visually represented the distribution of specific genes within gene modules in distinct subpopulations, accompanied by annotations outlining the functions of these gene modules. D Subsequently, macrophages were classified into M1 and M2 subtypes based on the marker of M1 subtype macrophages, denoted as m1up1 for high marker expression and m1dn2 for low marker expression in M1 subtype cells Role of Aggrephagy related genes-mediated endothelial in GB Tumor-associated endothelial cells (TAEC) have been recognized as pivotal players in GB pathogenesis and progression. They facilitate critical pathways, notably angiogenesis, and contribute to the recruitment of immune cells from the bloodstream to the tumor tissue [[109]35, [110]36]. In this investigation, we meticulously isolated 2647 TAECs from TC within the scRNA-seq data. Our objective was to elucidate the influence of Aggrephagy on the function of endothelial cells in the context of GB. To explore the intricate relationship between Aggrephagy and endothelial cells, we identified and characterized seven distinct endothelial NMF clusters that demonstrated positive expression of specific ARGs: DYNC1H1+TAEC-C1, TUBA1C+TAEC-C2, PARK7+TAEC-C3, TUBB4B+TAEC-C4, TUBA1A+TAEC-C5, UBB+TAEC-C6, and Unclear-TAEC-C7 (Fig. [111]5A). Our pseudotime trajectory analysis unveiled a progressive up-regulation of genes such as TUBA1A, DYNC1H1, PARK7, DYNLL1, and UBB during Aggrephagy-mediated cell development. Conversely, genes including TUBA4A, UBC, DYNLL2, and TUBA1C exhibited significant down-regulation over time (Fig. [112]5B). Furthermore, an examination of intercellular communication between these clusters unveiled the fibroblast cluster mediated by DYNC1H1, displaying the highest intensity in cellular communication pathways (Fig. [113]5C, D). This highlights the potential significant role of DYNC1H1-expressing endothelial cells in facilitating cellular interactions. Figure [114]5E portrays the outgoing and incoming signaling patterns between the ARGs cell clusters and macrophages. Notably, signaling molecules like SPP1, VISFATIN, and ANGPT were prominently activated. In the SCENIC analysis, we observed substantial heterogeneity in the expression activity of transcription factors (TFs) within the seven NMF clusters (Fig. [115]4F). The DYNC1H1+TAEC-C1 group exhibited the highest TF activity, while TUBA1C+TAEC-C2 and TUBB4B+TAEC-C3 retained higher activity compared to the remaining four groups. Furthermore, TFs such as JUN_extended, FOSB_extended, MEF2A_extended, and SOX4_extended demonstrated greater variability in activity between the DYNC1H1+TAEC group and the UBB+TAEC group (Fig. [116]5F). Finally, Fig. [117]5G illustrated the differential average expression of common signaling pathway genes in the ARGs cluster. Significantly, the UBB+TAEC-C6 and TUBB4B+TAEC-C4 groups exhibited heightened expression of pro-inflammatory factors. Fig. 5. [118]Fig. 5 [119]Open in a new tab Analysis of ARGs-regulated intercellular communication in endothelial cells, incorporating pseudotime and SCENIC analyses, produced several key insights. A Endothelial cell ARGs subtypes were visually represented using the UMAP method. B Pseudotime analysis of endothelial cells unveiled the dynamics of ARGs expression in TAEC. C The total number of intercellular communications between endothelial cells and macrophages was quantified. D The communication weights between endothelial cells and macrophages, as well as between endothelial cells themselves, were visually presented. E A depiction of signal output versus signal input for ARGs-regulated endothelial cell communication, with color depth indicating communication intensity. F SCENIC analysis yielded a heatmap showing the activity of transcription factors. G Another heatmap illustrated the expression of genes associated with immunosuppression Ucell scoring analysis for Aggrephagy related genes-mediated endothelial Furthermore, we conducted pathway enrichment analysis specifically for tumor-associated endothelial cells (TAEC). In Fig. [120]6A, we employed three distinct methods for enrichment analysis of Aggrephagy-related TAEC: Ucell, singscore, and ssgsea. The outcomes revealed that, across all three enrichment methods, the DYNC1H1-positive clusters and the TUBA1C-positive clusters exhibited notable up-regulation in functional activity. Consequently, we chose to proceed with the Ucell method for further analysis. Subsequently, we scrutinized the shared differential genes among these seven cell clusters. The results unveiled that TUBB4B-positive clusters and PARK7-positive clusters had the least overlap with other clusters. In contrast, the DYNC1H1 clusters displayed the most extensive intersection with the other clusters, indicating their heightened functional activity (Fig. [121]6B). We then visualized these findings through heatmaps with pathway annotations, illustrating that both DYNC1H1-positive clusters and TUBA1C-positive clusters exhibited significant functional up-regulation (Fig. [122]6C). Notably, key pathways such as EMT, TGF-beta pathway, and the unfolded protein response were prominently enriched. Our pathway activity was projected onto the UMAP map of endothelial cells categorized by ARGs. A comparison with Fig. [123]6A highlighted that all of these pathways were concentrated within two types of endothelial cells, namely TUBA1C+ and DYNC1H1+ (Fig. [124]6D). Simultaneously, we observed that UBB-positive clusters, as a whole, displayed marked downregulation in functional activity (Fig. [125]6E). In Fig. [126]6F, we compared the expression of key pathways in DYNC1H1+ and TUBA1C+ clusters. The results unveiled that DYNC1H1+ clusters were significantly enriched in the angiogenesis pathway. Moreover, this cell class exhibited substantial upregulation in pathways indirectly linked to angiogenesis, such as the TGF-beta signaling pathway, hedgehog signaling pathway, glycolysis, and mitotic spindle formation. While TUBA1C-positive clusters and DYNC1H1-positive clusters shared certain pathways, including glycolysis and MOTRC1, the DYNC1H1+ isoforms displayed heightened activity compared to the TUBA1C+ isoforms in numerous critical pathways, such as the TGF-beta signaling pathway and EMT, among others. Fig. 6. [127]Fig. 6 [128]Open in a new tab Functional enrichment analysis of endothelial cells under ARG regulation uncovered critical insights. A We employed three distinct methods to evaluate the overall functional activity of each endothelial cell subtype. Ucell was chosen for subsequent analysis. B Venn analysis unveiled overlapping functional enrichments among endothelial cell subtypes governed by ARGs. C A heatmap depicted the results of functional enrichment analysis among these subtypes, with the number of asterisks indicating the significance level. In the “Direction” column, blue denoted down-regulation, while red indicated up-regulation. D Notably, essential pro-tumorigenic pathways, including aggrephagy-related pathways and angiogenesis, were predominantly distributed in the DYNC1H1+TAEC and TUBA1C+TAEC cell subpopulations, as projected on the UMAP visualization. E A density distribution plot demonstrated the prevalence of pathways promoting tumor cell metastatic proliferation. F Differences in pathways between TUBA1C+TAEC and DYNC1H1+TAEC subtypes were visualized using half-violin plots, featuring violin plots on the left and box plots on the right Identify center cells for cell communication and clarify the intact cell-communication chains Upon the identification of Aggrephagy-associated endothelial cells and macrophages, our endeavor was to delve into the intricate web of cellular communication between each specific cell subtype. Our objective was to pinpoint the key subtypes that would unveil the comprehensive communication network intertwined with the central subtype. Within the realm of cellular communication between endothelial cells and macrophages, we embarked on a meticulous exploration, dissecting their roles as signaling sources. What transpired was a compelling revelation—among endothelial cells, clusters exemplified by DYNC1H1+ and TUBA1C+ exhibited extensive and prolific interactions with their cellular counterparts. Simultaneously, among macrophages, the HSP90AA1+ and TUBA1C+ clusters emerged as pivotal entities, fostering extensive interactions with other cells (Fig. [129]7A, B). The visualization of pathway-related heatmaps unfurled an interesting dimension. DYNC1H1+ endothelial cells and HSP90AA1+ macrophages shone with significantly high scores in pathways of relevance. Notably, DYNC1H1+ endothelial cells made a direct contribution to pathways associated with cancer (Fig. [130]7C). Consequently, we pivoted our focus toward DYNC1H1+ endothelial cells as the epicenter of our investigation into the intricate network of cellular communication pathways. Within this framework, we defined cells that transmitted signals to DYNC1H1+ endothelial cells as upstream cells, and when DYNC1H1+ endothelial cells emitted signals, these signals were deemed ligands, with the cells receiving these signals designated as downstream cells. This conceptual framework unveiled clusters of ligand-releasing cells upstream of endothelial cells and clusters of receptor-expressing cells downstream. Notably, in the realm of interactions involving DYNC1H1+ endothelial cells, PARK7+TAEC, TUBA1C+TAEC, and DYNC1H1+TAEC emerged as the most active in releasing ligands among upstream cells. As for downstream cells expressing receptors, they exhibited marked activity, with the exception of the PARK7+TAEC subtype, hinting that TUBA1C+TAEC and the central cell itself acted as both ligands and receptors in the communication network (Fig. [131]7D). A foray into the KEGG pathways linked with these Ligand-Receptor (LR) pairs unearthed intriguing findings. Among the most active pathways, those tethered to tumors such as Focal adhesion, PI3K-Akt pathway, ECM-receptor interaction, and gap junction came to the fore. Cells associated with these pathways predominantly belonged to the DYNC1H1+TAEC and TUBA1C+TAEC subpopulations (Fig. [132]7E). Figure [133]7F offered a panoramic view of the complete communication chain involving central cells. Seven pivotal up-regulated pathways took center stage, and an assortment of ligand-receptor pairs contributed to these pathways. The results underscored the significant role of receptors such as FLT1, PDGFRB, KDR, and ANGPT2 in multiple pathways. The pathways intertwined with endothelial cells revealed a complexity and diversity that surpassed their macrophage counterparts. Fig. 7. [134]Fig. 7 [135]Open in a new tab We have identified central hub communication cells within ARGs-regulated endothelial and macrophage populations, shedding light on their upstream and downstream communication chains. A The communication results using ARGs-regulated endothelial cells as signal outputs are visualized via hierarchical diagrams. B A similar visualization is provided for cellular communication using ARGs-regulated macrophages as signal outputs. C The heatmap showcases the strength of communication pathways for each cell subtype, highlighting the top 10 up-regulated pathways for each cell population. DYNC1H1 subtype endothelial cells, exhibiting the highest activity, were selected as central hub cells to delineate the communication chain. D In DYNC1H1+TAEC cell communication, we pinpointed the receptors of the central cell receiving signals from upstream cells and the ligands secreted by the central cell transmitting signals to downstream cells. Color and dot size represent p values and fold changes. E The most active communication pathways were examined, encompassing both upstream and downstream networks. F We meticulously outlined the complete upstream and downstream communication chain of DYNC1H1+TAEC, highlighting the pivotal role of aggrephagy pathways intertwined with significant tumor-promoting pathways. In the pathway annotation, color and length of the bars signify the mean pathway score differences and t values, while dot color and size represent −log10(P) values and log2FC, respectively. The central column “Pathway” elucidates the central role played by the DYNC1H1+ subtype in seven key pathways, along with its comprehensive communication chain involving specific ligand receptors Role of monocytes regulated by ARGs in TME We meticulously curated a cohort of 656 monocytes nestled within tumor clusters and subsequently subjected them to an NMF downclassing process based on ARGs, which, in turn, unveiled a trifecta of monocyte subtypes: TUBA1A+, TUBA1B+, and the Unclear cohort (Fig. [136]8A). Intriguingly, a comprehensive cellular communication analysis ensued between these identified subtypes and macrophages, unearthing a stark contrast. The TUBA1B+ subtype monocytes asserted dominance, surpassing their counterparts in both the number of communication instances and the overall communication potency (Fig. [137]8B, C). Delving further into the intricate web of interactions between these diverse isoforms and macrophages, a heatmap visualization of signal output and input was conducted. This endeavor unveiled disparities among the three subpopulations, most notably in key signaling pathways such as GALECTIN and TGF-beta. However, the monocyte population, as a collective entity, exhibited a surplus of input signals compared to output signals, with the TUBA1B+ subpopulation emerging as the most significant contributor (Fig. [138]8D). Subsequently, we embarked on an in-depth exploration of metabolic pathway enrichment within these three subtypes. The findings underscored the TUBA1B+ subtype’s marked up-regulation in metabolic pathways, encompassing sphingolipid metabolism, pyruvate metabolism, glycolipid metabolism, and fatty acid metabolism (Fig. [139]8E). The mosaic was completed with an examination of the disparities in the expression of key genes within the TME. The outcome illuminated the high expression of TUBA1B+ isoforms in factors intrinsically linked to pathways such as Extracellular Matrix (ECM), Matrix Metalloproteinases (MMPs), and Transforming Growth Factor-beta (Fig. [140]8F). Fig. 8. [141]Fig. 8 [142]Open in a new tab Analysis of monocytes under the regulation of ARGs reveals crucial insights. A ARGs-mediated monocytes are vividly depicted using the UMAP method. B We gauge the significance of cellular communication by assessing the weights of interactions between ARGs-regulated monocytes and macrophages. C We tally the number of these vital cellular exchanges between ARGs-regulated monocytes and macrophages. D We provide an illustrative overview of the ligand receptors engaged in input and output signals, facilitating the communication between ARGs-regulated monocytes and macrophages. E We delve into metabolic enrichment, shedding light on the distinct metabolic characteristics of ARGs-regulated monocytes. F We scrutinize the expression patterns of pivotal tumor-promoting genes within monocytes, revealing their regulation by ARGs ARGs-mediated TME patterns predict prognosis in Bulk RNA-seq Merging the TCGA matrix data from GB with clinical information, we embarked on the endeavor of aggrephagy scoring in this bulk dataset. Initially, we scrutinized the global expression disparities within the aggrephagy gene set between tumor and normal samples, revealing an intriguing revelation: the overall scoring in normal samples surpassed that in their tumor counterparts (Fig. [143]9A). Subsequently, we meticulously cataloged the distinctive marker genes associated with each cellular subpopulation linked to ARGs in the single-cell data. With this knowledge in hand, we proceeded to perform single-cell scoring predicated on the expression of these genes. Our focus remained on discerning the scoring distinctions between normal and tumor samples for each distinct cellular subtype. The outcome underscored that in most cell subpopulations, tumor samples exhibited higher scores than their normal counterparts. Particularly noteworthy were endothelial cells, including DYNC1H1+, TUBA1C+, TUBA1A+, PARK7+, and UBB+, which exhibited the most significant disparities. Macrophages, such as DYNLL1+, TUBA1C+, and PARK7+, were also noteworthy in this regard. Notably, monocytes displayed no significant score differences between normal and tumor samples, while HSP90AA1+ macrophages demonstrated higher scores in normal samples compared to their tumor counterparts (Fig. [144]9B). A subsequent pivotal step involved conducting Kaplan–Meier survival analysis (Fig. [145]9C). The results unearthed that certain cellular subtypes, such as TUBA1A+ monocytes, TUBA1B+ monocytes, TUBA1C+ TAEC, UBB+ TAEC, HSP90AA1+ macrophages, DYNLL1+ macrophages, TUBA1C+ macrophages, TUBB4B+ macrophages, and UBB+ macrophages, exhibited a grave prognosis. To fortify our findings, we extended our investigation to GB data within the GEO, consolidating and deconvoluting it before repeating the K–M survival analysis. This endeavor yielded significant results, indicating that solely TUBA1B+ monocytes portrayed an unfavorable prognosis. In the realm of endothelial cells, TUBA1A+ and PARK7+ cells exhibited a favorable prognosis, while DYNC1H1+, TUBA1C+, UBB+, and TUBB4B+ cells endured a dire prognosis. As for macrophages, HSP90AA1+, TUBA1C+, PARK7+, TUBB4B+ exhibited a bleak prognosis (Fig. [146]10A). Our quest for deeper insights led us to engage in logistic regression analysis, culminating in the creation of informative dot plots. These revelations pointed to the stability in prognosis across two datasets for specific cell subtypes, including PARK7+ macrophages, TUBA1B+ monocytes, TUBA1C+ macrophages, TUBA1C+ TAEC, and TUBB4B+ macrophages, as well as UBB+ TAEC. Significantly, all of these subtypes exhibited an unfavorable prognosis (Fig. [147]10B). Fig. 9. [148]Fig. 9 [149]Open in a new tab An in-depth prognosis of ARGs within the TCGA dataset unveils valuable insights. A We examine the differential expression of ARGs in normal and tumor samples within the TCGA-GB database, revealing distinctive patterns showcased through boxplots. B We conduct an evaluation of differential scoring for cell subtypes influenced by ARGs. This scoring method is applied to TCGA samples, considering the expression of marker genes linked to ARGs-regulated cellular subtypes in bulk data. We compare the distinctions between normal and tumor samples, with the number of asterisks denoting statistical significance. C Our findings from Kaplan–Meier survival analysis illustrate survival probabilities, where the yellow curve represents low-scoring ARG subtypes, while the red curve signifies high-scoring ARG subtypes Fig. 10. [150]Fig. 10 [151]Open in a new tab Prognostic evaluation of ARGs in GEO data, along with their associated cellular subtypes. A Depicting the outcomes through Kaplan–Meier survival analysis, vividly represented by K–M curves. B A comparative assessment of the prognostic impact of individual ARG subtypes in TCGA and GEO datasets. The color of the data points signifies the magnitude of risk, quantified as logHR, while the size of the data points indicates the level of significance, expressed as −log10P ARGs-mediated TME patterns affects immunotherapy in GB We delved further into the influence of ARGs-associated cellular subtypes on GB immunotherapy. Within the TCGA dataset, DYNC1H1+TAEC and TUBA1A+TAEC exhibited higher scores in the immune-unresponsive group compared to the immune-responsive group. Conversely, in the GEO database, subtypes TUBA1C+TAEC, TUBA1A+TAEC, TUBA1A+mono, TUBA1C+mac, and PARK7+mac demonstrated elevated scores in the immune-unresponsive group in contrast to the immune-responsive group (Fig. [152]11A, B). A thorough logistic regression analysis yielded logOR values less than 0 for subtypes DYNLL1+mac, PARK7+mac, TUBA1A+TAEC, TUBA1C+mac, and TUBA1C+TAEC in both the TCGA and GEO datasets, signifying a potential risk in the context of immunotherapy (Fig. [153]11C). Fig. 11. [154]Fig. 11 [155]Open in a new tab Immunotherapy prognosis regarding ARGs-regulated cell subtypes. A Evaluation of scores for ARGs-regulated cell subtypes in TCGA data, distinguishing non-responsive from responsive groups. B Assessment of scores for ARGs-regulated cell subtypes in GEO data, distinguishing non-responsive from responsive groups. C Identification of risk or protective factors in immunotherapy through logistic regression analysis. A lower logOR value indicates a risk factor, while a higher logOR value indicates a protective factor. The size of the data points reflects their significance, measured as −log10P Molecular docking analysis for TUBA1C and UBB In our investigation, we deliberately selected ARGs associated with cell subtypes of paramount importance in prognosis and immunotherapy, focusing on their sensitivity to various chemotherapeutic agents. Notably, TUBA1C exhibited pronounced sensitivity to a spectrum of chemotherapeutic agents (Fig. [156]12A). Given the prominence of TUBA1C and UBB expression in tumor samples within the TCGA dataset, we embarked on an extensive analysis of these two candidate proteins. To illuminate their potential roles in drug therapy, we conducted a meticulous molecular docking analysis. Employing Autodock Vina v.1.2.2, we computed binding energies for each interaction. The results are compelling, revealing that the drug candidates engage both protein targets through conspicuous hydrogen bonding and robust electrostatic interactions. Significantly, the two candidate proteins, TUBA1C and UBB, manifest low binding energies of −8.862 and −7.789 kcal/mol, indicative of highly stable binding (Fig. [157]12B, C). Fig. 12. [158]Fig. 12 [159]Open in a new tab Drug sensitivity analysis and molecular docking analysis. A Drug sensitivity analysis for key ARGs expressed in TME cells. B Molecular docking analysis for TUBA1C interacting with Quercetin. C Molecular docking analysis for UBB interacting with Quercetin Discussion As the intricate world of selective autophagic mechanisms, such as aggrephagy, unfolds, there’s a growing interest in unraveling their connection to cancer. Nevertheless, in the realm of gliomas, comprehensive investigations into the role of aggrephagy within the TME remain scarce. This groundbreaking study embarks on an exploration of the single-cell aggrephagy landscape within this context, indicating its value for both prognosis and immunotherapy. In meticulous detail, we scrutinized the functional enrichment, metabolic intricacies, and more, as governed by ARGs in individual cells. Furthermore, we delved into the intricate network of information exchange between the mesenchymal stromal cells in the TME and the immune cells. Our findings underscore the pivotal significance of these distinctive cellular subpopulations in shaping both prognosis and the landscape of immunotherapy, as validated through an analysis of bulk RNA-seq data. In the context of Homo sapiens, phagocytes within the brain typically comprise two components: tissue-resident microglia and bone marrow-derived macrophages (BMDMs) [[160]37]. As described by Sevenich et al., the onset of gliomas triggers a notable transformation in the composition of phagocytes within the tumor core. This shift entails a transition from a prevailing presence of microglia to an abundance of TAMs, while the peritumoral tissue remains akin to normal tissue, predominantly characterized by microglial presence [[161]38]. In our study, we similarly observed this conspicuous alteration in cell proportions preceding and following tumorigenesis. Existing studies have indicated that an escalation in the proportion of TAMs portends a grim prognosis. TAMs have been demonstrated to directly facilitate tumor metastasis and indirectly promote neo-angiogenesis. In high-grade GB, TAMs exhibit immunosuppressive attributes [[162]9, [163]39–[164]41]. Simultaneously, we noted a significant increase in both the quantity and proportion of endothelial cells among the mesenchymal stromal cells in the tumor group. Notably, endothelial cells are known to undergo proliferation and differentiation during angiogenesis, resulting in an elevated count and heightened functional activation. They play a pivotal role in tumor development [[165]42]. Hence, aligning our findings with the existing literature, our focus gravitated toward macrophages and endothelial cells in the TME. We embarked on a journey to pinpoint aggrephagy subtypes based on gene expression, embarking on an in-depth exploration of the regulatory influence of ARGs on cellular metabolism. Our quest was to unearth the distinctive roles played by these cells within the TME. Additionally, we ventured further to probe the intricate interactions between each subtype and the array of cells coexisting in the TME. By weaving together the regulatory networks inherent to each cell’s functions and their cross-communication, we unraveled the profound significance of these isoforms in shaping both prognosis and the landscape of immunotherapy in GB. In our initial macrophage analysis, we conducted a comparative study with microglia, pinpointing differentially expressed genes. Among the down-regulated genes in macrophages, several hold pivotal roles as tumor suppressors. TSC22, for instance, acts as a tumor suppressor by inhibiting the TGF-beta signaling pathway, which can otherwise promote epithelial-mesenchymal transition [[166]43, [167]44]. Another significant gene, CACNA1A, serves as a subunit of voltage-dependent calcium channels and has demonstrated importance in tumor immune checkpoint blockade therapy, exerting a notable tumor-suppressive effect in conditions like lung cancer [[168]45, [169]46]. Additionally, PDK4, marked by its down-regulation, is recognized as a poor prognostic indicator in cancer, boasting crucial anti-tumor effects [[170]47]. Conversely, among the genes upregulated in macrophages, a subset assumes roles as key oncogenes. Notably, members of the S100A family, such as S100A10, S100A9, S100A8, and S100A6, exhibit substantial upregulation across various cancer types, indicating their oncogenic potential [[171]48, [172]49]. Studies have demonstrated that S100A10 can stimulate malignant tumor development by activating the mTOR-Signaling Pathway [[173]50]. Moreover, research by Sarvaiya et al. has unveiled the role of CXCL1, CXCL2, and CXCL3 in promoting tumor growth in cancers like pancreatic cancer and melanoma [[174]51]. Notably, TAMs in gliomas have been found to exhibit significant CXCL3 enrichment in multiple studies [[175]52, [176]53]. In light of these collective findings and our own analysis, we contend that the macrophages characterized in this study represent a substantial risk factor in gliomas. Further investigation into the regulatory mechanisms of ARGs in the TME, particularly concerning TAMs, is warranted. In the realm of macrophages, subsequent to identifying distinct NMF subtypes, we uncovered noteworthy variations in metabolic activity. The Unclear subpopulation and the HSP90AA1 subpopulation exhibited significantly heightened metabolism. The Unclear subtype, specifically characterized by the upregulation of one or more aggrephagy genes (adjusted p value <0.05, log2FC > 0 but <1), emerged as the most metabolically active. Notably, aggrephagy-associated subtypes exhibited active metabolic profiles, albeit generally weaker than the Unclear subtype. This underscores the bidirectional influence of ARGs on macrophage metabolism, ultimately favoring upregulation in metabolic pathways. In vitro experiments have suggested that TAMs may support tumor cells by synthesizing fatty acids, among other functions. In more malignant tumors, pathways such as fatty acid synthesis, oxidation, and glycolysis in TAMs have demonstrated significant upregulation. Glycolysis, in particular, leads to lactic acid production, promoting the secretion of Vegf and arg1 by TAMs, consequently driving tumor differentiation [[177]54]. Considering these findings, we postulate that cell subtypes displaying substantial metabolic activity, including HSPAA190+, TUBA1C+, PARK7+, and the Unclear subtype, play pivotal roles in fueling tumor growth and metastasis. Notably, HSP90AA1 has been identified as an immunosuppressive factor in cancer [[178]55]. Our exploration of macrophage functionally enriched pathways unveiled that HSP90AA1+ subtypes display pronounced activity in pathways related to the response of misfolded proteins. This analysis, coupled with metabolic activity findings, suggests that HSP90AA1 may regulate the expression of HSP90AA1 in ARGs to maintain TAM homeostasis, aiding in the phagocytosis or repair of misfolded proteins during their immunosuppressive functions. However, a notable discovery was that macrophages of the HSP90AA1+ subtype exhibited significant downregulation of the chemokine-mediated signaling pathway, particularly the CXCL3 gene. CXCL3, known to be expressed in TAMs, is recognized for promoting the recruitment of TAMs and subsequent infiltration into the TME. Simultaneously, CXCL3 promotes neutrophil infiltration and angiogenesis [[179]56, [180]57]. Intriguingly, the CCL3 gene within the HSP90AA1+ subtype showed significant upregulation, and in the TME, CCL3 plays a role in mediating the selective recruitment of TAMs and M2 macrophages [[181]58, [182]59]. Integrating these findings, it becomes evident that while TAMs downregulate non-selective TAMs chemokines, they concurrently upregulate selective TAMs chemokines, ultimately promoting TAM infiltration into the TME. In stark contrast, the TUBA1C+ subtype showed significant upregulation in the CXCL3 pathway and insignificant downregulation in the CCL3-associated pathway, underscoring its role as a significant promoter of tumorigenesis and development. The UBB+ subtype exhibited the second-highest upregulatory activity in the chemokine-mediated signaling pathway, following closely behind the TUBA1C+ subtype. This subtype also displayed substantial upregulation in genes like B4GALT1 and TLR2. B4GALT1, recognized as a poor prognostic marker in gliomas [[183]60], likely contributes to its immunosuppressive properties against the tumor [[184]61]. On the contrary, TLR2 plays a role in regulating TAM function and promoting tumorigenesis [[185]62]. Furthermore, protumor genes, including HSP90AA1, TUBA1C, UBB, PARK7, and TUBB4B, are prominently expressed as marker genes in ARGs [[186]63, [187]64], suggesting that ARGs collectively act to promote tumorigenesis in macrophages. We proceeded to scrutinize the expression of ARGs in both M1 and M2 macrophage subtypes, revealing a substantial increase in ARGs expression within the M1 subtype compared to non-aggrephagy counterparts. Conversely, no significant difference was observed in the M2 subtype. Previous research has delineated the M2 subtype as a promoter of immunosuppressive factors, while the M1 subtype typically champions antitumor factors [[188]65]. Our comprehensive analysis suggests that ARGs play a pivotal role in fostering TAM infiltration and dampening immunity within macrophage regulation. Consequently, they exhibit heightened expression in M1 macrophages, thereby inhibiting the antitumor activity of M1 subtypes and indirectly augmenting the proportion of immunoinactive TAMs that lack the immune response in the TME. We similarly delineated aggrephagy subtypes within endothelial cells. Through comprehensive enrichment analysis heatmaps and scatter plots depicting pathway activity distribution, we unveiled the prominence of functionally active subtypes like DYNC1H1+ and TUBA1C+ in key pathways pivotal to tumorigenesis and development. The DYNC1H1+ subtype exhibited direct upregulation of the angiogenesis pathway, closely linked with hypoxia, glycolysis, and other critical routes that contribute to angiogenesis [[189]66]. Furthermore, the DYNC1H1+ subtype indirectly stimulated angiogenesis by elevating these pathways. According to Porter et al., the glycolytic pathway in endothelial cells significantly fosters tumor proliferation and anti-apoptosis [[190]67]. Additionally, the hedgehog pathway emerged as a substantial pro-tumorigenic factor, especially in brain cancer, where sonic hedgehog is a recognized drug therapy target [[191]68]. Notably, this subtype was a direct promoter of epithelial-mesenchymal transition, a pathway profoundly implicated in determining GB aggressiveness [[192]69]. The TGF-beta pathway, indirectly encouraging EMT in GB, was also notably upregulated [[193]70]. These insights underscore the DYNC1H1+ subtype as a conspicuous tumor-promoting factor within the glioma TME. Conversely, the TUBA1C+ subtype exhibited significant upregulation of pro-tumorigenic pathways like fatty acid synthesis and glycolysis. However, it also demonstrated heightened activity in pathways such as p53 and TNFA, directly suppressing tumor metastasis. On the contrary, other subtypes associated with aggregated autophagy displayed a general downregulation of pathway activities, including p53, TNFA, and other anti-tumor pathways. This collectively suggests their contribution as risk factors in GB. Notably, all other subtypes exhibited considerably less activity than the DYNC1H1+ subtype in the key pathways of angiogenesis and EMT, emphasizing the pivotal role of the DYNC1H1+ subtype as the foremost pro-tumorigenic factor within endothelial cells. Moreover, an even more intriguing revelation manifests as the direct upregulation of aggrephagy-related pathways in the DYNC1H1+ subtype. To begin with, the unfolded protein signaling pathway serves as a guardian, shielding endothelial cells from the clearance of misfolded proteins. Furthermore, the PI3K-AKT-MTOR pathway, a quintessential aggrephagy-related pathway, not only fosters the formation of protein aggregates but also triggers the autophagy-related pathway for the efficient elimination of aggregated proteins, thereby preserving the delicate homeostasis of endothelial cells. It is noteworthy that the MTORC1 pathway, situated as a downstream branch of the PI3K-AKT-MTOR pathway, has been substantiated to exhibit activity in IDHWT GB cells [[194]71]. This compelling revelation also unravels that the KRAS pathway, positioned upstream of the PI3K-AKT pathway, indirectly augments the aggrephagy-related pathway. This pivotal insight suggests that ARGs transcend their role in merely regulating endothelial cells to induce risk factors within the TME. Moreover, ARGs extend a protective mantle over the most detrimental DYNC1H1 subtype by orchestrating the aggrephagy-related pathway within these at-risk subtypes. This orchestration enables at-risk subtypes to effectively eliminate erroneous proteins during pivotal processes like angiogenesis, EMT, and various pathways while maintaining cellular homeostasis. Consequently, this robust coordination perpetuates a tumor-promoting role within the TME. Notably, we also observed analogous protective mechanisms within the TUBA1C+ subtype. The DYNC1H1+ subtype not only exhibited pronounced activity in functional enrichment pathways but also displayed heightened intercellular communication, particularly with macrophages, surpassing all other subtypes. Consequently, we designated it as the central cell for analysis, aiming to unravel the pathways through which this cell subtype communicates within the TME. Our exploration extended to identifying the upstream and downstream cells involved in these pathways, along with the ligand receptors governing this communication network. The top 20 pathways, closely linked to tumorigenesis, were meticulously selected for investigation. Upon scrutinizing the pathway heatmap, we observed that the DYNC1H1+ subtype scored significantly in the “pathways in cancer.” Notably, it exhibited robust involvement in pathways that promote glioma tumorigenesis. In particular, the DYNC1H1+ subtype of endothelial cells demonstrated extensive upstream cell communication, implying its ability to enhance functionality through autocrine regulation. Thus, it was positioned as the central hub cell to explore the interconnected cellular communication network. Within this network, we noted that it exhibited significant connections with other endothelial cell subtypes in pathways like gap junction and focal adhesion. Gap junctions are known to facilitate cancer cell adhesion, migration, and carcinogenesis. Moreover, they play a vital role in angiogenesis, notably driven by VEGF in the nervous system. In our study, we deduced that the upregulation of the gap junction pathway by DYNC1H1+ endothelial subtypes was primarily aimed at establishing effective gap junction intercellular communication (GJIC), thereby fostering an efficient interaction network. GJIC establishment is a key promoter of endothelial cell proliferation. Furthermore, we observed that platelet-derived growth factor receptor (PDGFR) served as the receptor associated with the gap junction pathway. PDGF, in turn, induced the expression of VEGF in tumor endothelial cells, promoting their proliferation, differentiation, and the recruitment of pericytes, consequently inducing neoangiogenesis. The focal adhesion pathway, activated predominantly by VEGF induction, plays a crucial role in maintaining endothelial cell adhesion during angiogenesis by enhancing the PI3K-Akt pathway. In addition to these endothelial cell-specific signaling pathways, we discovered the significance of TGF-beta signaling in macrophages. This pathway, initiated by communication from central hub cells, plays a pivotal role in promoting EMT and immunosuppression in TAMs. Remarkably, central hub cells exhibited profound engagement within the PI3K-Akt pathway, forging connections with various endothelial cells and macrophage subtypes via angiogenesis-promoting receptors like ANGPT2, EDN1, and PDGFB. This implies that DYNC1H1+ endothelial cells, acting as central players, foster extensive functional interactions with neighboring cells. Their concerted efforts in bolstering aggrephagy-related pathways serve to collectively uphold homeostasis, thereby perpetuating a stable progression of tumorigenesis within the TME. We also conducted NMF subtyping of monocytes, revealing that nearly all subtypes exhibited robust cellular interactions with macrophages. However, a distinct subtype, TUBA1B+, stood out with the highest metabolic activity among monocytes. It displayed significant enrichment in critical pathways, including ECM and neo-angiogenesis, and exhibited strong communication with macrophages through pathways such as TGF-beta and MIF. This underscores the role of the TUBA1B+ subtype in promoting neoangiogenesis within myeloid cells, particularly by fostering the execution of pro-tumorigenic processes in the TME. These findings shed light on the dynamic interactions among myeloid cell subtypes, specifically highlighting the influence of the TUBA1B+ subtype through hypoxia-related signaling pathways. To validate our scRNA-seq findings, we incorporated bulk data into our analysis. Notably, in the differential analysis of ARGs, the normal group exhibited significantly higher expression than the tumor group. However, when we examined the scores of ARGs-regulated mesenchymal stromal cells, a different trend emerged. In most subtypes, the tumor group displayed notably higher scores compared to the normal group. Only monocytes and certain macrophage subtypes showed no significant differences. This intriguing observation implies that the regulatory influence of ARGs on cells within the TME plays a pivotal role in glioma development. The impact of ARGs extends to the metabolism, functionality, and intercellular communication of critical mesenchymal cells and immune cells within the glioma microenvironment. Endothelial cells and macrophages, once under ARG regulation following tumorigenesis, become activated and collaboratively enhance angiogenesis, tumor migration, and immunosuppression pathways. Cox univariate analysis further underlines the risk factors presented by almost all cell subtypes regulated by ARGs in the context of GB. Kaplan–Meier curves illustrate that these cell subtypes are associated with poor prognosis in both TCGA and GEO datasets. Notably, cell subtypes with adverse outcomes, such as TUBA1B+mono-C2, TUBA1C+TAEC-C2, TUBB4B+TAEC-C4, UBB+TAEC-C6, TUBA1C+Mac-C3, PARK7+Mac-C4, and TUBB4B+Mac-C5, exhibit extensive intercellular communication with the central cell, DYNC1H1+TAEC. In summary, the influence of ARGs on NMF cell subpopulations is linked to a significant poor prognosis in GB, underscoring their vital role in shaping the survival outcomes of GB patients. We proceeded to substantiate the impact of ARGs-mediated cell subtypes on ICB immunotherapy. Notably, our analysis revealed that TUBA1A+TAEC subtypes exhibited significant immunosuppression in both TCGA and GEO datasets. Concurrently, TUBA1C-mediated cells displayed immunosuppressive characteristics in both TAMs and TAECs, which are pivotal mesenchymal cell types within the glioma microenvironment. These findings underscore the substantial role played by the microtubule protein family in GB. As supported by Zhu et al., TUBA1C emerges as a key prognostic factor in gliomas, with elevated expression in GB relative to LGG and normal tissues [[195]72]. This consistency is further reinforced by immunohistochemical analyses in the study by Hu et al. Consequently, TUBA1C holds promise as a potential prognostic marker and therapeutic target. The detrimental role of TUBA1C in both TCGA and GEO datasets is unequivocal, whereas the function of the UBB gene in GB remains enigmatic. Several studies have reported its overexpression in other carcinomas, including gastric cancer, pancreatic ductal adenocarcinoma, and cervical cancer. In these investigations, Scarpa et al. observed that UBB exerts a substantial pro-survival effect on tumor cells in gastric cancer, participating in processes such as hyperproliferation, migration, and invasion. Yang et al. indicated a close relationship between UBB and base repair excision, while Peng et al. found that down-regulating UBB significantly decreases ubiquitin levels, essential for cancer cell growth [[196]73–[197]76]. By analyzing the docking of the 2 proteins with a certain chemotherapeutic agent, we infer that the 2 proteins can act as targets for chemotherapy in GB. In summary, the conclusions of these studies consistently support UBB’s role in enhancing the stability of the intracellular environment of tumor cells and promoting their migration and invasiveness. It is essential to acknowledge the preliminary nature of our study, which entails certain limitations. Firstly, our scRNA-seq data originates from one GEO dataset, constraining both the sample size and sequencing depth. To enhance the robustness of our findings, expansion to a larger patient cohort is imperative. Furthermore, we observed that several pathways upregulated by the DYNC1H1+TAEC subtype, such as TNFA-SIGNALING and the P53 pathway, exhibit tumor-suppressive roles. Additionally, the APOPTOSIS pathway exerts negative regulation in angiogenesis concerning endothelial cells. Consequently, further exploration is warranted to elucidate the delicate balance between tumor-promoting and tumor-suppressing factors, as well as to unravel the underlying mechanisms that drive ARGs toward a more pronounced tumor-promoting role. Additionally, the limited presence of control group samples in the TCGA-GB database may introduce potential errors into the analytical outcomes. The integration of actual tissue samples would bolster the validation of our predictive results. Conclusion In our investigation, we delved into the distinctive roles and intricate interaction networks of ARGs-regulated cellular subtypes within the TME using scRNA-seq data. Our exploration unearthed pivotal subtypes displaying a notable upregulation of aggrephagy-related pathways. Subsequently, we substantiated the impact of each subtype on GB prognosis and immunotherapy through bulk RNA-seq analysis. In essence, our findings underscore the capacity of ARGs to activate specific TME cells, driving processes like immunosuppression, immune evasion, and tumor neoangiogenesis, ultimately fostering tumor progression and aggravating patient outcomes. Moreover, we observed that ARGs also function to safeguard particular high-risk cell subtypes, thus exacerbating poor prognosis among GB patients by preserving TME cellular equilibrium. Acknowledgements