ScRNA-seq data from 44 primary glioma samples was obtained from a previous study (18). To ensure data quality, single cells expressing fewer than 500 expressed genes, over 20% mitochondrial transcripts, or containing more than 50% ribosomal transcripts were excluded from further analysis. Additionally, we removed genes that were expressed in fewer than three single cells. We used the DoubletFinder Python package was utilized to identify potential doublets (19). The filtered dataset consisted of 22,8156 cells, which were analyzed using (20). We normalized gene expression using Seurat’s LogNormalize method, applying a scale factor of 10,000 (21). Highly variable genes (n=2000) were selected and their expression values were scaled prior to principal component analysis (PCA). Batch effects were corrected using the Harmony R package (22). Data analysis was performed using functions from the Harmony and Seurat R packages, including NormalizeData, FindVariableFeatures, ScaleData, RunPCA, FindNeighbors, FindClusters, and RunUMAP. Cell cycle phases were scored using Seurat’s CellCycleScoring function (23).

To assess the score of large-scale chromosomal copy number variations (CNVs) in somatic cells, the inferCNV R package (24) was utilized. A raw counts matrix, an annotation file, and a gene/chromosome position file were prepared according to the specified data prerequisites (https://github.com/broadinstitute/inferCNV). T/NK cells and B cells were subsequently designated as the reference cells for the analysis.

We developed a novel framework that encompasses five distinct algorithms—CCAT (16), CytoTRACE (17), Monocle3 (25), PAGA (26), and Slingshot (27)—each meticulously designed to assess differentiation capacity of cells using sscRNA-seq data. CCAT, based on the concept of entropy-rate, was executed through the SCENT R package (16). In contrast, CytoTRACE provides an unsupervised framework for predicting relative differentiation states from single-cell transcriptomes, utilizing the CytoTRACE R package (17). An increase in scores calculated by both CCAT and CytoTRACE indicates a higher degree of cellular differentiation. Pseudotime analysis was conducted using the Monocle3 R package, the Slingshot R package, and the PAGA method in the SCANPY Python package (27). The framework code can be accessed at the following GitHub repository: https://github.com/Caolab2024/Cancer_stem_cells/tree/main.

Cell type abundances were estimated from the bulk kidney expression data in the TCGA cohort using the BisqueRNA R package (28). A PCA-based method was employed to deconvolute the seven primary kidney cell types based on scRNA data for subsequent analyses.

We applied SCENIC (29), a novel computational approach for inferring regulatory networks and identifying TFs from scRNA-seq data, to individual cells. Subsequently, we employed receiver operating characteristic curve (ROC) analysis to identify regulons that were preferentially expressed in distinct cell clusters based on transcription factors or their target genes.

We obtained spatial transcriptome data from 28 specimens available in the OSF repository(https://osf.io/4q32e/), comprising a total of 88,793 spots. Using the Seurat, an R package for single-cell genomics (20), we performed the following steps: (1) normalized and scaled the UMI counts using SCTransform and identified the most variable features; (2) reduced dimensionality and clustered the spots with RunPCA, applying the default parameters and the top 30 principal components.

To deconvolve the transcriptome data into cell-type-specific gene expression profiles, we utilized RCTD (30), a computational method for resolving cell types in complex biological samples. RCTD can help uncover the cellular composition, function and interactions in biological research.

We estimated cell-type dependencies using MISTy (31), a method for inferring mutual information between cell types. MISTy was applied to the RCTD estimates from all slides, and utilizing a multi-view model with a parameter that weighted the estimates of neighboring cell types (effective radius = 15 spots). We interpreted the median standardized importances of each view across all slides as indicators of cell-type colocalization or mutual exclusion in various spatial contexts.

CellChat analysis was performed using the CellChat R package (32). This tool enables the inference and analysis of cell-cell communication networks from single-cell RNA sequencing data. By utilizing CellChat, we were able to explore intercellular signaling pathways and identify key communication patterns among distinct cell populations.

We searched several public databases, including The Cancer Genome Atlas (TCGA, http://portal.gdc.cancer.gov/),ArrayExpress(https://www.ebi.ac.uk/biostudies/arrayexpress), Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn/), and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) for datasets that met the following criteria: (1) more than 60 samples per cohort; (2) Affymetrix Human Genome U133 Plus 2.0 Array or high-throughput sequencing platforms; and (3) primary tumors from patients who did not receive any treatments before resection. We obtained a total of 2285 samples from 6 cohorts: TCGA-GBMLGG (n = 683) (33), CCGA1 (n = 413) (34), CCGA2 (n = 273) (34), GSE16011 (n = 262) (35), GSE108474 (n = 490) (36), and E-MTAB-3892 (n =164) (37). Transcripts per million (TPM) data for TCGA-GBMLGG were downloaded from UCSC Xena database. For data generated using the Affymetrix Human Genome U133 Plus 2.0 Array, we applied the robust multiarray averaging (RMA) algorithm from the Affy R package for preprocessing. Finally, We log2 transformed, z-score normalized, and removed batch effects from the gene expression data across all cohorts using the surrogate variable analysis (SVA) algorithm (38).

We sought to develop a precise and robust GSCS for predicting glioma patient prognosis. To achieve this, we built an artificial intelligence network incorporating 429 algorithm combinations, integrating 27 algorithms from traditional regression, machine learning, and deep learning approaches. These algorithms included CoxTime, DeepSurv, DeepHit, Logistic-Hazard, PC-hazard, Akritas, Coxboost, VSOLassoBag, RSF, GBM, SuperPC, obliqueRSF, CForest, GLMBoost, BlackBoost, Rpart, Survreg, Ranger, Ctree, LASSO, plsRcox, survival-SVM, Ridge, Enet, XGBoost, Boruta, and stepwise Cox. Initially, we conducted univariate Cox regression to identify prognostic GSC markers in the TCGA cohort at a significance level of P < 0.05. Subsequently, we applied the 429 algorithm combinations to these markers to develop predictive models within the TCGA cohort. The predictive performance of each algorithm combination was evaluated using C-indices across all validation cohorts. The optimal algorithm combination was selected based on the highest mean C-index. We stratified glioma patients into high- and low-risk groups according to the optimal cutoff value obtained by the survminer R package. The network code can be accessed at the following GitHub repository: https://github.com/Caolab2024/Cancer_stem_cells/tree/main.

We conducted gene set variation analysis (GSVA) and gene set enrichment analysis (GSEA) using the MSigDB database employing the GSVA and clusterprofiler (39, 40) R packages. The differentially expressed genes (log2FC > 1, adjusted P <0.05) between the low- and high-risk groups were subjected to Metascape for enrichment analysis (41).

Cell transfection was conducted using two distinct siRNAs targeting TUBA1C, synthesized by Ribobio (Guangzhou, China) and delivered with Lipofectamine 3000 (Invitrogen, USA). The siRNA sequences are provided in Supplementary Table 1 . RNA extraction from tissues or cell lines was performed using TRIzol (Thermo) followed by cDNA synthesis with the PrimeScript™RT kit preceded gene expression quantification was then carried out via SYBR qPCR Master Mix on the Roche LightCycler 480 (Roche, GER) (42). Primer sequences were obtained from Tsingke Biotech (Beijing, China) and detailed in Supplementary Table 1 .

In each well of the 96-well plates, two thousand treated cells were seeded, followed by the addition of the CCK-8 labeling reagent for further processing (43). Observations and assessments were conducted on days 0, 1, 2, 3, 4, and 5 to monitor the cell responses and outcomes.

Transfected LN299 and U87 cells (siNC, siTUBA1C-1, siTUBA1C-2 groups) were seeded at a density of 1,000 cells per well in a 6-well plate and incubated for 14 days to allow colony formation. The medium was replaced every 3 days. After incubation, the media were aspirated, and the cells were washed with phosphate-buffered saline (PBS). The cells were then fixed with 4% paraformaldehyde (PFA) for 20 minutes at room temperature, followed by staining with 0.5% crystal violet (Solarbio, China) for 20 minutes (44). After staining, excess crystal violet was rinsed off with water, and once dry, the number of colonies per well, defined as clusters of at least 50 cells, was counted. Colony numbers were compared across the siNC, siTUBA1C-1, and siTUBA1C-2 groups to evaluate the impact of gene silencing on colony formation.

Following the transfection process, once cell confluence reached 95%, the transfected cells were seeded into 6-well plates. A sterile 200 μL pipette tip was used to draw a straight line, facilitating the gentle removal of unattached cells and debris with PBS. Subsequently, serum-free cell medium was replenished to sustain the cell culture. Photographs were captured at both the 0-hour and 48-hour time points at identical locations for comparative analysis.

Cells were seeded at a density of 2×10⁴ per well in 200 μL of serum-free medium within the upper chamber, which was either coated or left uncoated with matrix glue from BD Biosciences, USA. The lower chamber contained 700 μL of 10% complete medium. After a growth period of 36 hours, the cells were fixed, stained, and photographed for quantification.

All mice were housed in the animal facility of Northern Jiangsu People’s Hospital Affiliated to Yangzhou University and maintained in a Specific-Pathogen Free (SPF) environment. The animal experiments were approved by the Ethics Committee of Northern Jiangsu People’s Hospital Affiliated to Yangzhou University. LN299 NC and Si-1 cell lines were cultured to the logarithmic growth phase, washed twice with PBS, and collected. The cell concentration was adjusted to 1×10⁷ cells/ml. A 100 μl cell suspension was subcutaneously injected into the right flank of each 6-week-old female nude mouse (BALB/c-nu), using 5 mice per group. Tumor length and width were measured every three days with calipers, and body weight was recorded. Tumor volume was calculated using the formula: Volume (mm³) = 0.5 × Length (mm) × Width (mm)². On day 21 post-injection, the mice were sacrificed, and the tumors were excised and weighed.

The excised tumor tissues were fixed in 10% neutral formalin for 24 hours, followed by paraffin embedding and sectioning at a thickness of 4 μm. The sections were deparaffinized in xylene, rehydrated through a graded ethanol series, and endogenous peroxidase activity was blocked with 3% hydrogen peroxide. Antigen retrieval was performed using citrate buffer (pH 6.0) in a pressure cooker. After cooling to room temperature, the sections were blocked with goat serum for 30 minutes and then incubated overnight at 4°C with a Ki67 primary antibody (1:200 dilution, Abcam). The following day, the sections were washed three times with PBS for 5 minutes each, incubated with a secondary antibody for 30 minutes, developed with DAB, counterstained with hematoxylin, dehydrated, and mounted.

Data analysis, statistics, and plotting were performed using R 4.3.1. Continuous variables were assessed using the Wilcoxon rank-sum test or the T test. The optimal cut-off value was determined with the survminer R package. Survival analysis was conducted using Cox regression and Kaplan-Meier methods via the survival R package. The pROC package was employed to implement the ROC curve for predicting binary categorical variables. The timeROC R package calculated the time-dependent area under the curve (AUC) for survival variables. Unless otherwise specified, P < 0.05 is considered statistically significant.

The study flow diagram is shown in Figure 1 . We analyzed 22,8156 cells from 44 samples that passed quality control steps to uncover the cellular and molecular heterogeneity of cancer cells in human gliomas. Unsupervised clustering identified 38 clusters with distinct gene expression patterns ( Figures 2A–D ). InferCNV analysis was employed to distinguish malignant and non-malignant cells based on the CNV score ( Figure 2E ). Each cluster was assigned to a cell type using InferCNV analysis and marker gene expression ( Figures 2F, G ).

After conducting additional dimensionality reduction and clustering of glioma cells, our bespoke analytical framework enabled us to pinpoint GSCs, as indicated in Figure 3A . Utilizing the PAGA algorithm, we demonstrated the differential potential of GSCs as they evolve into various tumor cell subpopulations ( Figure 3B ). Analyses integrating Monocle3 and the Slingshot algorithm further elucidated the fundamental role that GSCs subpopulations play as the origin of tumoral evolution ( Figures 3C, D ). Both the CCAT and CytoTRACE algorithms highlighted the pronounced stem-like qualities inherent within GSCs ( Figures 3E, F ). Then, we performed a comprehensive enrichment analysis of the highly expressed marker genes across various tumor subpopulations, utilizing the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) databases. Our findings revealed significant enrichment of GSC marker genes within pathways associated with cell proliferation, particularly those involving the cell cycle ( Figure 3G ). To further elucidate the proliferative capacity of GSCs, we conducted a detailed assessment of their cell cycle profiles, which demonstrated a predominant presence in the G2M and S phases ( Figure 3H ). This observation underscores the pivotal role of proliferation in the maintenance of stemness (45). Moreover, our analysis uncovered a striking prevalence of GSCs within samples of higher malignancy, including recurrent GBM, and gliomas with wild-type (WT) IDH ( Figure 3H ).

Utilizing bulk RNA-seq data from the TCGA cohort, we evaluated the prognostic significance of five distinct tumor cell populations. Initially, we identified cell-type-specific marker genes at the single-cell level, characterized by a log2 fold change exceeding 1 and an adjusted P below 0.05, and proceeded to determine their hazard ratios (HRs) based on overall survival (OS), progression-free interval (PFI), disease-specific survival (DSS), and disease free interval (DFI) within the TCGA cohort. Our observations revealed that GSC marker genes exhibited the most elevated HRs for OS, PFI, DFI, and DSS ( Figures 4A–D ). We ascertained that these GSCs constitute independent prognostic indicators for OS, PFI, and DSS (P < 0.05) ( Figures 4E–G ). Nevertheless, GSCs did not emerge as independent risk factors for DFI, a finding potentially attributable to the paucity of DFI data ( Figure 4H ). Kaplan-Meier analysis further demonstrated that a higher abundance of GSCs was associated with adverse outcomes across OS, PFI, DFI, and DSS (P < 0.05) ( Figures 4I–L ). Moreover, the area under the curve (AUC) values for the prediction of OS, PFI, DSS, and DFI consistently surpassed 0.7 ( Figures 4I–L ), underscoring the formidable prognostic predictive power of GSCs.

As previous studies (46) showed, TFs are key for cell fate specification. Each cell type exhibited a unique TFs pattern ( Figure 5A ). Strong enrichment was in E2F1, E2F2, E2F7, and BRCA1 regulon activity in GSCs ( Figure 5B ). Figure 5C displays the expression patterns of these four TFs, which correlated with worse prognosis in glioma patients from the TCGA cohort (P < 0.05) ( Figure 5D ). These findings indicated that the four TFs might be essential for maintaining the stemness of GSCs.

We reduced the dimensionality and clustered the spatial transcriptome data, identified the predominant cell types at each spot based on the back-convolutional via the scRNA-seq reference data, and chose the two samples with the highest GSCs abundance for the next step in the analysis ( Figures 6A, B ). Based on MISTy results, GSCs had a close spatial location to myeloid-derived suppressor cells (MDSCs) ( Figure 6C ). Then, we found that the MIF pathway was active between GSCs and MDSCs (Figure 6D). This result confirmed the previous conclusion that GSCs activate MDSCs to suppress immune responses by secreting MIF (47).

To further quantify the abundance of GSCs using key genes and improve the ability to predict prognosis in gliomas, we developed a GSCS on a novel artificial intelligence network. We first performed univariate Cox analysis to select the most specific GSCs marker (log2 fold change > 1, adjusted P = 0) with prognostic value (P < 0.05). We then fitted 429 algorithm combinations on the TCGA cohort and computed the C-index for each combination on the validation cohorts. The combination of VSOLassoBag and RSF had the highest mean C-index of 0.764 ( Figure 7A ). VSOLassoBag identified 26 genes, which were used by RSF to construct the GSCS ( Supplementary Figures 1A, B ). We stratified glioma patients into high- and low-risk groups based on the optimal cutoff from the TCGA cohort. Our results showed that the high-risk group had significantly worse OS than the low-risk group in all cohorts (P < 0.05) ( Figures 7B–G ). Moreover, time-dependent ROC curves demonstrated the robust and stable performance of the GSCS in the all cohorts ( Figures 7B–G ).

The GSCS outperformed age, grade, gender, IDH status, MGMT promoter status, karnofsky score (KPS), 1p/19q co-deletion,TP53 protein expression, and recurrent status in terms of the C-index across the all cohorts ( Figure 8A ). We also compared the GSCS with other pulished signatures in the TCGA, GSE108474, GSE16011, E-MTAB-3892, CGCA1 and CGCA2 cohorts ( Figures 8B–H ). The GSCS had the highest C-index among all signatures in the all cohorts.

We calculated the GSCS score for 33 cancers from the TCGA ( Figure 9A ) and divided them into High-risk and Low-risk groups. We observed that high GSCS scores correlated with worse prognosis across all 15 cancers (P < 0.05) ( Figure 9B ). These results suggest that GSCS may also be a key factor affecting the prognosis of various other cancers.

To explore the biological mechanisms underlying the association between GSCS and proliferative features, we performed pathway analysis on the GSCS score. We found that the score is strongly correlated with several tumorigenic pathways, such as the G2M DNA replication checkpoint, angiogenesis, E2F targets, and epithelial-mesenchymal transition (EMT) (P < 0.05) ( Figure 10A ). We also observed significant differences in proliferation-related pathways between the two risk groups (P < 0.05) ( Figure 10B ). The DEGs between the low- and high-risk groups were enriched in immune-related and proliferation-related pathways (P < 0.05) ( Figure 10C ). Furthermore, GSEA of kyoto encyclopedia of genes and genomes (KEGG) terms revealed that the high-risk group was enriched for ecm-receptor interaction, cell cycle, P53 signaling pathway, and DNA replication (P < 0.05) ( Figure 10D ). These findings also demonstrated the critical role of proliferation in maintaining stemness (48).

To investigate the key role of TUBA1C in the pathogenesis of glioma, we selected the LN299 and U87 glioma cell lines as research subjects and used targeted siRNA-mediated knockdown technology to regulate the expression of TUBA1C ( Figure 11A ). We evaluated the impact of TUBA1C knockdown on the viability of glioma cells using the CCK-8 assay. The experimental results showed that the OD value of glioma cells decreased significantly after TUBA1C expression was silenced ( Figures 11B, C ). To further investigate the impact of TUBA1C suppression on the migration and invasion ability of glioma cells, we conducted in-depth analyses. Using the Transwell migration and invasion assay, we observed that the invasion and migration ability of glioma cells were significantly reduced after TUBA1C knockdown ( Figures 11D–F ). Clonogenic assays confirmed that the proliferation ability of glioma cells was weakened after TUBA1C knockdown ( Figures 11G, H ). The wound healing experiment provided quantitative data, showing that TUBA1C knockdown significantly slowed the wound healing process, indicating that cell migration ability was significantly inhibited ( Figures 11I, J ). Taken together, we have revealed the key role of TUBA1C in promoting the proliferation, migration, and invasion of glioma cells. This discovery underscores the importance of TUBA1C as a potential therapeutic target in glioma treatment and may provide new ideas and methods for future glioma treatment.

We injected LN299 cell lines transfected with either TUBA1C-NC or TUBA1C-Si-1 subcutaneously into nude mice and dynamically observed the body weight and tumor growth of the mice. Through qPCR analysis, we compared the differences in TUBA1C expression between the two groups and confirmed the effectiveness of gene silencing in the xenograft tumors ( Supplementary Figure 2 ). The results indicated that, compared to the NC group, the tumor weight and volume were significantly reduced in the Si-1 group mice, while there was no significant difference in body weight between the groups ( Figures 12A–D ). Immunohistochemical analysis of the mouse tissues showed that the percentage of Ki67 positive cells in the tumor tissues of the Si-1 group mice was also significantly reduced, suggesting that Si-1 treatment effectively inhibits tumor growth and cell proliferation ( Figures 12E, F ). In summary, these results highlight the potential of TUBA1C as a therapeutic target in glioma treatment strategies.