Lei Cao1
Xu Lu
Hao Wu
Xiaye Miao- 1Department of Oncology, The Affiliated Suqian First People’s Hospital of Nanjing Medical University, Suqian, China
- 2Department of Oncology, The Affiliated Huai’an Hospital of Xuzhou Medical University and the Second People’s Hospital of Huai’an, Huai'an, China
- 3Department of Laboratory Medicine, Northern Jiang su People's Hospital, Yangzhou, Jiangsu, China
Background: Gliomas are aggressive brain tumors associated with a poor prognosis. Cancer stem cells (CSCs) play a significant role in tumor recurrence and resistance to therapy. This study aimed to identify and characterize glioma stem cells (GSCs), analyze their interactions with various cell types, and develop a prognostic signature.
Methods: Single-cell RNA sequencing data from 44 primary glioma samples were analyzed to identify GSC populations. Spatial transcriptomics and gene regulatory network analyses were performed to investigate GSC localization and transcription factor activity. CellChat analysis was conducted to infer cell-cell communication patterns. A GSC signature (GSCS) was developed using machine learning algorithms applied to bulk RNA sequencing data from multiple cohorts. In vitro and in vivo experiments were conducted to validate the role of TUBA1C, a key gene within the signature.
Results: A distinct GSC population was identified, characterized by high proliferative potential and an enrichment of E2F1, E2F2, E2F7, and BRCA1 regulons. GSCs exhibited spatial proximity to myeloid-derived suppressor cells (MDSCs). CellChat analysis revealed an active MIF signaling pathway between GSCs and MDSCs. A 26-gene GSCS demonstrated superior performance compared to existing prognostic models. Knockdown of TUBA1C significantly inhibited glioma cell migration, and invasion in vitro, and reduced tumor growth in vivo.
Conclusion: This study offers a comprehensive characterization of GSCs and their interactions with MDSCs, while presenting a robust GSCS. The findings offer new insights into glioma biology and identify potential therapeutic targets, particularly TUBA1C, aimed at improving patient outcomes.
Introduction
Gliomas are the most common and aggressive primary tumors affecting the central nervous system. According to data from the Central Brain Tumor Registry of the United States (CBTRUS) for the years 2013 to 2017, these neoplasms accounted for approximately 25% of all adult primary brain tumors and 81% of malignant central nervous system tumors in the United States (1). The National Comprehensive Cancer Network (NCCN) Guidelines classify glioma as a diverse group of neoplasms. Ranging from low-grade gliomas (LGGs), such as surgically treatable pilocytic astrocytomas to highly invasive and virtually incurable glioblastoma multiforme (GBM) (2). Despite extensive research into molecular therapies targeting oncogenic pathways and immune checkpoints in gliomas, significant improvements in patient outcomes have remained elusive. This situation underscores the necessity for continued investigation into novel therapeutic approaches for this challenging group of tumors (3, 4).
Glioma stem cells (GSCs), characterized by their stem cell attributes, constitute a minor subset within the larger glioma cell population. The majority of the glioma mass comprises differentiated progeny, a characteristic conferred upon GSCs due to the capability for self-renewal (5). Within the cancer stem cells (CSCs) milieu, immune cells such as cancer-associated fibroblasts (CAFs), tumor-associated macrophages (TAMs), and myeloid-derived suppressor cells (MDSCs) secrete cytokines like TGFβ significantly contributing toward the EMT-mediated invasion of CSCs (5, 6). Contemporary research indicates a significant role of CSCs in glioma recurrence and resistance to chemoradiotherapy (6, 7). Targeting CSCs and their supportive microenvironment has emerged as a promising strategy for developing novel and more effective treatments for gliomas, with the aim of improving patient outcomes and addressing the challenge of recurrence. Flow cytometry or magnetic techniques facilitate the enrichment of CSC populations from bulk tumors, leveraging specific cell surface markers for selection. However, the expression of specific markers is not consistently observed across all glioma stem cells. It is essential to recognize that CSCs, when exposure to novel microenvironments, are susceptible to alterations in their state, a phenomenon observed regardless of the methodology used for CSC enrichment. Notably, in vitro culture conditions can prompt variations in surface marker expression and modulate the intrinsic biological states of glioma cells (8). Consequently, examining CSCs without prior selection provides a valuable opportunity to investigate their inherent characteristics, potentially yielding insights into the processes by which CSCs originate and differentiate into the cells implicated in gliomagenesis.
The enhanced resolution for cell-type identification and characterization, previously unattainable with bulk RNA sequencing (RNA-seq), has been facilitated by Single-cell RNA sequencing (scRNA-seq) (9, 10). Numerous researchers have already leveraged the advantages of single-cell sequencing to successfully discover cancer biomarkers and identify potential therapeutic targets (10–14). A primary objective of scRNA-Seq research is to decipher the hierarchical differentiation structures of complex tissues (15). To achieve this, it necessitates an unbiased measure of the differentiation potential of individual cells, thus enabling the recognition of either stem or multipotent progenitor cells and the establishment of a ranking of individual cells along potency gradients of differentiation (16). However, current in silico methods present challenges in differentiating between adult stem cells with long-term regenerative capabilities and more differentiated cells. Although models based on gene expression possess the potential to surmount these constraints (17), the extent of their applicability across varied developmental systems and an array of single-cell sequencing methodologies remains to be fully elucidated.
In this investigation, we developed a novel framework consisting of five distinct algorithms specifically designed to identify GSCs using scRNA-seq data. Furthermore, we utilized the SCENIC algorithm and spatial transcriptomics (ST) analysis to elucidate the transcription factor (TF) activities and cell communications within these GSCs. Bulk RNA-seq deconvolution highlighted the significant role of GSCs in predicting poor patient prognosis. Based on a comprehensive combination of 429 algorithmic combinations, we established a GSC Signature (GSCS). Ultimately, we conducted in vivo and in vitro experiments to empirically validate the malignant characteristics associated with TUBA1C, a gene that identified as the most critical component within the GSCS.
Methods
Collection and pre-processing of scRNA-seq data
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).
InferCNV analysis
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.
Identification of GSCs
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.
Abundance of cell types in bulk RNA-seq data
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.
Gene regulatory network analysis
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.
Processing of glioma spatial transcriptome sequencing data and inferring cellular localization
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.
Collection and pre-processing of bulk RNA-seq data
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).
Development of signatures using an artificial intelligence network
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.
Functional annotation of the GSCS
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 and RT-qPCR
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.
Cell counting
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.
Colony formation
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.
Wound healing
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.
Transwell
Cells were seeded at a density of 2×104 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.
Subcutaneous tumor xenograft in nude mice
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×107 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.
Ki67 immunohistochemistry staining
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.
Statistical analysis
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.
Results
Clustering and cell-type identification of single-cell RNA-seq data
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).
Figure 1. An illustration of the general workflow adopted in this study.
Figure 2. Clustering and cell-type identification of single-cell RNA-seq data. (A–F) UMAP projections of 228,156 aggregate single cells from 44 primary glioma samples showing the composition of different cell types in human gliomas. UMAP projections are shown by cluster numbers (A), by the patient (B), by pathological type (C), by grade (D), by Copy Number Variation (CNV) score (E) and by cell types (F). (G) Dot plot showing marker gene expression and cell cycle score for different cell types.
Identification and characterization of glioma stem cells
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).
Figure 3. Identification and characterization of glioma stem cells (GSCs). (A) UMAP projections of 92,014 aggregate glioma cells are shown by cell annotation. (B–D). Trajectory inference using PAGA (B), Slingshot (C), and Monocle3 (D). (E, F) UMAP projections of 92,014 aggregate glioma cells are shown by Cytotrace score (E) and CCAT score (F). (G) Marker genes and Kyoto Encyclopedia of Genes and Genomes (KEGG)/Gene Ontology (GO) enrichment analysis for each glioma cell subpopulations. (H) Differential distribution of glioma cell subsets across cell cycle phases (left), histopathological classifications (middle), and IDH mutation status (right).
GSCs correlated with unfavorable prognosis
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.
Figure 4. Higher GSCs abundance in glioma, linked with poor prognosis. (A–D) GSCs abundance’s association with poor overall survival (OS) (A), disease-specific survival (DSS) (B), progression-free interval (PFI) (C), and disease free interval (DFI) (D) manifested through HR values for cell type marker genes from Cox proportional hazards regression. Each dot is a gene, with HR value on x-axis and cell type on y-axis. (E–H). Multivariate Cox analyses to identify the risk factors in OS (E), DSS (F), PFI (G), DFI (H). (I–L). Time-dependent receiver operating characteristic (ROC) analysis (top) and Kaplan-Meier (KM) curves (below) of the GSCs abundance in OS (I), DSS (J), PFI (K), DFI (L).
Unique TF profile associated with 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.
Figure 5. Unique transcription factor (TF) activity associated with GSCs. (A) Heatmap showing differences in TF activity scored by SCENIC. (B) TF activity of E2F1, E2F2, E2F7 and BRCA1 projected on UMAP. (C) TF expression of E2F1, E2F2, E2F7 and BRCA1 projected on UMAP. (D) Kaplan–Meier survival curves for E2F1, E2F2, E2F7 and BRCA1.
Spatial transcriptomics of glioma reveals cellular localization
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).
Figure 6. Spatial transcriptomics of glioma reveals cellular localization and cell communication. (A) Spatial Transcriptomics-based UMAP clustering of sample 1 (top) and sample 2 (below). (B) Spatial localization of individual clusters in the sample 1 (top) and sample 2 (below). (C) The spatial distribution of different cell-types in the spatial transcriptomics reference calculated by RCTD in the sample 1 (top) and sample 2 (below). (D) MIF signaling pathway from CellChat results.
Construction of the GSCS
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).
Figure 7. An artificial intelligence network was utilized to develop and validate a consensus GSC signature (GSCS). (A) A total of 429 prediction models were developed using a 10-fold cross-validation framework, and the C-index of each model was computed across all datasets. (B–G). Kaplan-Meier curves and ROC curves of OS according to the GSCS in the (B) TCGA, (C) GSE16011, (D) GSE108474, (E) E-MTAB-3892, (F) CGG1, and (G) CGGA2.
Comparison of prognostic signatures
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.
Figure 8. (A) The C-index of the GSCS and other models developed in the (B) TCGA, (C) GSE16011, (D) GSE108474, (E) E-MTAB-3892, (F) CGCA1, and (G) CGCA2.
The GSCS exhibited a generalizability signature in pan-cancer
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.
Figure 9. Pan-cancer of GSCS. (A) GSCS score for 33 cancers in the TCGA database. (B) Kaplan-Meier curves of OS according to the GSCS in the TCGA database.
Potential biological peculiarities of the GSCS
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).
Figure 10. Biological peculiarities of the GSCS in the TCGA cohort. (A) MsigDB-based gene set variation analysis (GSVA) delineated the biological attributes of two risk groups. (B) t-Distributed Stochastic Neighbor Embedding (t-SNE) plots of kyoto encyclopedia of genes and genomes (KEGG) and reactome terms delineated the differences in pathway activity in the two risk groups. (C) Metascape-based enrichment analysis of differentially expressed genes between two risk groups. (D) Gene set enrichment analysis (GSEA) of KEGG terms for the GSCS. ***p < 0.001.
The silence of the TUBA1C inhibits the malignant biological behavior of glioma cells
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.
Figure 11. Silencing of TUBA1C Inhibits the Malignant Biological Behavior of Glioma Cells. (A) Schematic representation of the targeted siRNA-mediated knockdown of TUBA1C in LN299 and U87 glioma cell lines. (B, C) CCK-8 assay results showing a significant decrease in OD values of glioma cells following TUBA1C knockdown, indicating reduced cell viability. The results are presented as the mean ± SD of three independent experiments. (D–F) Transwell migration and invasion assays demonstrating a significant reduction in the migration and invasion abilities of glioma cells after TUBA1C knockdown. Quantitative analysis showed a significant decrease in the number of migrated and invaded cells. (G, H) Clonogenic assays confirming that the proliferation capacity of glioma cells is weakened upon TUBA1C knockdown. The number of colonies formed was significantly reduced compared to the control group. (I, J). Wound healing assays showing quantitative data that TUBA1C knockdown significantly slows the wound healing process, indicating a marked inhibition of cell migration capability. Wound closure percentage was significantly lower in TUBA1C knockdown cells. **p < 0.01; ***p < 0.001.
Silencing TUBA1C inhibits subcutaneous tumor growth in vivo
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.
Figure 12. Knockdown of TUBA1C Inhibits Subcutaneous Tumor Growth In Vivo. (A) Representative images of tumors from mice, showing the tumor size comparison between the NC group and the Si-1 group. Scale bar: 1 cm. (B) Comparison of tumor weights between NC and Si-1 groups of mice. (C) Body weight changes during the experiment, showing no significant difference between the two groups. (D) Tumor volume growth curves, showing that the tumor volume in Si-1 group mice was significantly smaller than that in the NC group. (E) Ki67 immunohistochemistry images showing Ki67 positive cells in tumors from NC and Si-1 groups. (F) Comparison of the percentage of Ki67 positive cells per field. ****p < 0.0001.
Discussion
Treatment for gliomas has remained unchanged since 2005, involving surgical resection, radiation and concurrent and adjuvant TMZ (49). Adjuvant TMZ marginally increased survival time for adults with gliomas; however, this agent has caused systemic toxicities and decreased the quality of life for patients (49, 50). Moreover, many tumors exhibit primary or acquired resistance to temozolomide (51). The anti-angiogenic antibody bevacizumab, which neutralizes VEGF, was approved for recurrent gliomas and increased PFS in patients (52). However, subsequent trials showed no improvement in OS for newly diagnosed patients (53). While bevacizumab normalized the tumor blood vessels and reduced symptoms such as oedema, it does not extend lifespan (54). Gliomas are widely infiltrative and resistant to standard therapies, remaining non-curative. Therefore, new therapeutic options are urgently needed. In this work, based on integrated analysis of scRNA, bulk RNA sequencing, ST and machine learning algorithms, we identified and characterized GSC and developed a GSCS based on maker genes from the GSC as a predictive model for glioma patient outcomes.
We identified a glioma cell type, GSCs, that demonstrate a predominant presence in the G2M and S phases of the cell cycle. Furthermore, we observed a strong enrichment of E2F1, E2F2, E2F7, and BRCA1 regulons in GSCs. These four TFs are implicated in proliferation and DNA repair, indicating the high proliferative potential of this subcluster. Previous studies have shown that these TFs regulate stemness (55, 56). Additionally, GSCs were enriched with genes previously linked to poor glioma outcomes, such as CENPF (57), TOP2A (58), NUSAP1 (59), PTTG1 (60), UBE2C (61), and UBE2S (62). Unlike previous studies that focused on individual glioma genes, we employed a single-cell decomposition approach to reveal the association of the GSC type proportion and reduced survival. Moreover, the GSCs encompass not only encompasses known glioma genes, but also novel targets, such as HMGN2, TUBB4B, and ARL6IP1, that warrant further investigation. In summary, our research confirms the role of previously known genes in glioma GSCs while also identifying novel targets and TF regulators that can enhance our understanding of this complex disease. With further investigation, these discoveries could pave the way for developing more effective therapeutic strategies for glioma patients.
We aimed to enhance our understanding of the factors affecting survival in glioma patients, given their poor prognosis. Our findings suggest that GSCs could serve as useful biomarker for guiding treatment and predicting outcomes. Current prognostic tools for glioma primarily rely on factors such as WHO grade, IDH mutations, 1p/19q codeletion, MGMT promoter methylation, TERT promoter mutations and EGFR amplification (63). Cell type markers can aid in comprehend cancer biology and supplement existing clinical practices. The marker genes associated with GSCs may pave the way for new expression-based prognostic technologies. As RNA sequencing technology has matured and clinical laboratories can now detect gene expression patterns with prognostic value (64). We developed an integrative pipeline to construct a GSCS using the expression profiles of GSCs marker genes. We validated this signature in six independent cohorts and confirmed that the best model is a combination of VSOLassoBag and RSF. The algorithm combinations eliminated low-value features, optimized the model, and enhanced its generalization ability (65). The GSCS performed well in predicting outcomes across all six cohorts, as shown by ROC and C-index analyses, indicating its potential clinical utility. To maximize the clinical utility of the GSCS, future research should focus on integrating this biomarker into routine diagnostic workflows. Investigating its application in personalized treatment plans could enhance therapeutic decision-making, particularly for patients with varying glioma subtypes. Additionally, exploring the GSCS in conjunction with existing prognostic factors may lead to more refined risk stratification models. Longitudinal studies assessing the signature’s predictive power in response to specific therapies will be crucial. Furthermore, expanding the signature’s validation across diverse populations and treatment settings will strengthen its relevance and applicability in clinical practice, ultimately improving patient outcomes.
Based on the RSF algorithm, TUBA1C was identified as the most significant gene in GSCS. TUBA1C is an isoform of α-tubulin that has been shown to play a critical role in the cell cycle and immune microenvironment of lung adenocarcinoma (LUAD). Elevated TUBA1C expression correlates with poor outcomes and with tumor-infiltrating immune cells (TIICs) in LUAD (66). Additionally, TUBA1C is upregulated in hepatocellular carcinoma (HCC) and pancreatic ductal adenocarcinoma (PDAC), where it predicts poor prognosis and enhances cell proliferation and migration (67, 68). Furthermore, a prior study indicated that TUBA1C was statistically associated with the expression of RP11-480I12.5 in breast cancer (BRCA) and demonstrated prognostic significance (69). TUBA1C has also been shown to promote aerobic glycolysis and cell growth via upregulation of YAP expression, thereby contributing to BRCA development. Our findings further revealed that TUBA1C knockdown significantly inhibited the malignant biological behaviors of glioma cells, demonstrating that TUBA1C is a promising target for the treatment of glioma.
This study has advanced our understanding of GSCs and their clinical relevance; however, we acknowledge several limitations. First, the cohorts exhibited had heterogeneity due to different in sequencing or microarray platforms. We harmonized the data using standard normal transformation, which was only partially effective. Second, we relied on retrospective samples, necessitating future validation in a prospective, large cohort. Third, we should conduct more in-depth and detailed molecular biology studies in both in vivo and in vitro experiments to uncover the molecular mechanisms of tumor recurrence and identify new therapeutic targets.
Conclusion
This study provides a comprehensive characterization of GSCs through integrated analysis of single-cell RNA sequencing, spatial transcriptomics, and machine learning approaches. We identified a distinct GSC population with high proliferative potential and developed a novel 26-gene GSCS that exhibits robust prognostic value across multiple cohorts. The signature demonstrated pan-cancer prognostic ability and an association with critical tumorigenic pathways. We validated the functional significance of TUBA1C, a key component of our signature, through in vitro and in vivo experiments. Silencing TUBA1C significantly inhibited glioma cell proliferation, migration, and invasion, as well as tumor growth in xenograft models. This study enhances our understanding of glioma biology and provides a clinically relevant prognostic tool and potential therapeutic targets.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material. The Raw data is available at: https://www.jianguoyun.com/p/DWH4qUEQ45m7Cxih8NwFIAA.
Ethics statement
The animal study was approved by The Ethics Committee Northern Jiangsu People’s Hospital Affiliated to Yangzhou University. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
LC: Funding acquisition, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. XL: Software, Visualization, Writing – review & editing. XW: Software, Validation, Writing – review & editing. HW: Conceptualization, Data curation, Methodology, Writing – original draft, Writing – review & editing. XM: Conceptualization, Data curation, Methodology, Writing – original draft, Writing – review & editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was supported by the Suqian Sci&Tech Program (Grant No. 32B0301).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1475235/full#supplementary-material
Supplementary Figure 1 | (A) Determination of the number of trees by minimizing error. (B) Variable importance of the top 26 genes determined using the random survival forest (RSF) algorithm.
Supplementary Figure 2 | Relative TUBA1C mRNA expression levels in xenograft tumors from si-NC and si-LN299 groups. qPCR analysis revealed that TUBA1C expression was significantly reduced in the si-LN299 group compared to the si-NC group. ***p < 0.001.
References
1. Ostrom QT, Patil N, Cioffi G, Waite K, Kruchko C, Barnholtz-Sloan JS. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2013-2017. Neuro Oncol. (2020) 22:iv1–iv96. doi: 10.1093/neuonc/noaa200
2. Nabors LB, Portnow J, Ahluwalia M, Baehring J, Brem H, Brem S, et al. Central nervous system cancers, version 3.2020, NCCN clinical practice guidelines in oncology. J Natl Compr Canc Netw. (2020) 18:1537–70. doi: 10.6004/jnccn.2020.0052
3. Wang X, Guo G, Guan H, Yu Y, Lu J, Yu J. Challenges and potential of PD-1/PD-L1 checkpoint blockade immunotherapy for glioblastoma. J Exp Clin Cancer Res. (2019) 38:87. doi: 10.1186/s13046-019-1085-3
4. Choma CT, Kaplan H. Bacillus thuringiensis crystal protein: effect of chemical modification of the cysteine and lysine residues. J Invertebr Pathol. (1992) 59:75–80. doi: 10.1016/0022-2011(92)90114-J
5. Suvà ML, Tirosh I. The glioma stem cell model in the era of single-cell genomics. Cancer Cell. (2020) 37:630–6. doi: 10.1016/j.ccell.2020.04.001
6. Prieto-Vila M, Takahashi RU, Usuba W, Kohama I, Ochiya T. Drug resistance driven by cancer stem cells and their niche. Int J Mol Sci. (2017) 18. doi: 10.3390/ijms18122574
7. Schulz A, Meyer F, Dubrovska A, Borgmann K. Cancer stem cells and radioresistance: DNA repair and beyond. Cancers (Basel). (2019) 11. doi: 10.3390/cancers11060862
8. Gisina A, Kholodenko I, Kim Y, Abakumov M, Lupatov A, Yarygin K. Glioma stem cells: novel data obtained by single-cell sequencing. Int J Mol Sci. (2022) 23. doi: 10.3390/ijms232214224
9. Ye B, Wang Q, Zhu X, Zeng L, Luo H, Xiong Y, et al. Chen T et al: Single-cell RNA sequencing identifies a novel proliferation cell type affecting clinical outcome of pancreatic ductal adenocarcinoma. Front Oncol. (2023) 13:1236435. doi: 10.3389/fonc.2023.1236435
10. Ye B, Hongting G, Zhuang W, Chen C, Yi S, Tang X, et al. Deciphering lung adenocarcinoma prognosis and immunotherapy response through an AI-driven stemness-related gene signature. J Cell Mol Med. (2024) 28:e18564. doi: 10.1111/jcmm.18564
11. Zhang P, Dong S, Sun W, Zhong W, Xiong J, Gong X, et al. Deciphering Treg cell roles in esophageal squamous cell carcinoma: a comprehensive prognostic and immunotherapeutic analysis. Front Mol Biosci. (2023) 10:1277530. doi: 10.3389/fmolb.2023.1277530
12. Zhang L, Cui Y, Mei J, Zhang Z, Zhang P. Exploring cellular diversity in lung adenocarcinoma epithelium: Advancing prognostic methods and immunotherapeutic strategies. Cell Prolif. (2024) 2024:e13703. doi: 10.1111/cpr.13703
13. Zhang P, Pei S, Zhou G, Zhang M, Zhang L, Zhang Z. Purine metabolism in lung adenocarcinoma: A single-cell analysis revealing prognostic and immunotherapeutic insights. J Cell Mol Med. (2024) 28:e18284. doi: 10.1111/jcmm.18284
14. Ye B, Ji H, Zhu M, Wang A, Tang J, Liang Y, et al. Single-cell sequencing reveals novel proliferative cell type: a key player in renal cell carcinoma prognosis and therapeutic response. Clin Exp Med. (2024) 24:167. doi: 10.1007/s10238-024-01424-x
15. Rozenblatt-Rosen O, Stubbington MJT, Regev A, Teichmann SA. The Human Cell Atlas: from vision to reality. Nature. (2017) 550:451–3. doi: 10.1038/550451a
16. Teschendorff AE, Maity AK, Hu X, Weiyan C, Lechner M. Ultra-fast scalable estimation of single-cell differentiation potency from scRNA-Seq data. Bioinformatics. (2021) 37:1528–34. doi: 10.1093/bioinformatics/btaa987
17. Gulati GS, Sikandar SS, Wesche DJ, Manjunath A, Bharadwaj A, Berger MJ, et al. Single-cell transcriptional diversity is a hallmark of developmental potential. Science. (2020) 367:405–11. doi: 10.1126/science.aax0249
18. Abdelfattah N, Kumar P, Wang C, Leu JS, Flynn WF, Gao R, et al. Single-cell analysis of human glioma and immune cells identifies S100A4 as an immunotherapy target. Nat Commun. (2022) 13:767. doi: 10.1038/s41467-022-28372-y
19. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. (2019) 8:329–337.e324. doi: 10.1016/j.cels.2019.03.003
20. Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. (2015) 33:495–502. doi: 10.1038/nbt.3192
21. Song G, Peng G, Zhang J, Song B, Yang J, Xie X, et al. Uncovering the potential role of oxidative stress in the development of periodontitis and establishing a stable diagnostic model via combining single-cell and machine learning analysis. Front Immunol. (2023) 14:1181467. doi: 10.3389/fimmu.2023.1181467
22. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. (2019) 16:1289–96. doi: 10.1038/s41592-019-0619-0
23. Tirosh I, Izar B, Prakadan SM, Wadsworth MH 2nd, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. (2016) 352:189–96. doi: 10.1126/science.aad0501
24. Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. (2014) 344:1396–401. doi: 10.1126/science.1254257
25. Cao J, Spielmann M, Qiu X, Huang X, Ibrahim DM, Hill AJ, et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature. (2019) 566:496–502. doi: 10.1038/s41586-019-0969-x
26. Wolf FA, Hamey FK, Plass M, Solana J, Dahlin JS, Göttgens B, et al. PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol. (2019) 20:59. doi: 10.1186/s13059-019-1663-x
27. Street K, Risso D, Fletcher RB, Das D, Ngai J, Yosef N, et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics. (2018) 19:477. doi: 10.1186/s12864-018-4772-0
28. Jew B, Alvarez M, Rahmani E, Miao Z, Ko A, Garske KM, et al. Accurate estimation of cell composition in bulk expression through robust integration of single-cell information. Nat Commun. (2020) 11:1971. doi: 10.1038/s41467-020-15816-6
29. Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. (2017) 14:1083–6. doi: 10.1038/nmeth.4463
30. Cable DM, Murray E, Zou LS, Goeva A, Macosko EZ, Chen F, et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nat Biotechnol. (2022) 40:517–26. doi: 10.1038/s41587-021-00830-w
31. Wu Y, Yang S, Ma J, Chen Z, Song G, Rao D, et al. Spatiotemporal immune landscape of colorectal cancer liver metastasis at single-Cell level. Cancer Discovery. (2022) 12:134–53. doi: 10.1158/2159-8290.CD-21-0316
32. Sun Z, Wang J, Zhang Q, Meng X, Ma Z, Niu J, et al. Coordinating single-cell and bulk RNA-seq in deciphering the intratumoral immune landscape and prognostic stratification of prostate cancer patients. Environ Toxicol. (2024) 39:657–68. doi: 10.1002/tox.23928
33. Brennan CW, Verhaak RG, McKenna A, Campos B, Noushmehr H, Salama SR, et al. The somatic genomic landscape of glioblastoma. Cell. (2013) 155:462–77. doi: 10.1016/j.cell.2013.09.034
34. Zhao Z, Zhang KN, Wang Q, Li G, Zeng F, Zhang Y, et al. Chinese glioma genome atlas (CGGA): A comprehensive resource with functional genomic data from chinese glioma patients. Genomics Proteomics Bioinf. (2021) 19:1–12. doi: 10.1016/j.gpb.2020.10.005
35. Gravendeel LA, Kouwenhoven MC, Gevaert O, de Rooi JJ, Stubbs AP, Duijm JE, et al. Intrinsic gene expression profiles of gliomas are a better predictor of survival than histology. Cancer Res. (2009) 69:9065–72. doi: 10.1158/0008-5472.CAN-09-2307
36. Gusev Y, Bhuvaneshwar K, Song L, Zenklusen JC, Fine H, Madhavan S. The REMBRANDT study, a large collection of genomic data from brain cancer patients. Sci Data. (2018) 5:180158. doi: 10.1038/sdata.2018.158
37. Kamoun A, Idbaih A, Dehais C, Elarouci N, Carpentier C, Letouzé E, et al. Integrated multi-omics analysis of oligodendroglial tumours identifies three subgroups of 1p/19q co-deleted gliomas. Nat Commun. (2016) 7:11263. doi: 10.1038/ncomms11263
38. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. (2012) 28:882–3. doi: 10.1093/bioinformatics/bts034
39. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118
40. Zhang S, Jiang C, Jiang L, Chen H, Huang J, Gao X, et al. Construction of a diagnostic model for hepatitis B-related hepatocellular carcinoma using machine learning and artificial neural networks and revealing the correlation by immunoassay. Tumour Virus Res. (2023) 16:200271. doi: 10.1016/j.tvr.2023.200271
41. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. (2019) 10:1523. doi: 10.1038/s41467-019-09234-6
42. Ma B, Qin L, Sun Z, Wang J, Tran LJ, Zhang J, et al. The single-cell evolution trajectory presented different hypoxia heterogeneity to reveal the carcinogenesis of genes in clear cell renal cell carcinoma: Based on multiple omics and real experimental verification. Environ Toxicol. (2024) 39:869–81. doi: 10.1002/tox.24009
43. Zhang H, Xia T, Xia Z, Zhou H, Li Z, Wang W, et al. KIF18A inactivates hepatic stellate cells and alleviates liver fibrosis through the TTC3/Akt/mTOR pathway. Cell Mol Life Sci. (2024) 81:96. doi: 10.1007/s00018-024-05114-5
44. Zhu C, Sun Z, Wang J, Meng X, Ma Z, Guo R, et al. Exploring oncogenes for renal clear cell carcinoma based on G protein-coupled receptor-associated genes. Discovery Oncol. (2023) 14:182. doi: 10.1007/s12672-023-00795-z
45. Kildisiute G, Kalyva M, Elmentaite R, van Dongen S, Thevanesan C, Piapi A, et al. Transcriptional signals of transformation in human cancer. Genome Med. (2024) 16:8. doi: 10.1186/s13073-023-01279-z
46. Chen K, Ma Y, Liu X, Zhong X, Long D, Tian X, et al. Single-cell RNA-seq reveals characteristics in tumor microenvironment of PDAC with MSI-H following neoadjuvant chemotherapy with anti-PD-1 therapy. Cancer Lett. (2023) 576:216421. doi: 10.1016/j.canlet.2023.216421
47. Otvos B, Silver DJ, Mulkearns-Hubert EE, Alvarado AG, Turaga SM, Sorensen MD, et al. Cancer stem cell-Secreted macrophage migration inhibitory factor stimulates myeloid derived suppressor cell function and facilitates glioblastoma immune evasion. Stem Cells. (2016) 34:2026–39. doi: 10.1002/stem.2393
48. Jiang H, Wei H, Wang H, Wang Z, Li J, Ou Y, et al. Zeb1-induced metabolic reprogramming of glycolysis is essential for macrophage polarization in breast cancer. Cell Death Dis. (2022) 13:206. doi: 10.1038/s41419-022-04632-z
49. Stupp R, Mason WP, van den Bent MJ, Weller M, Fisher B, Taphoorn MJ, et al. Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N Engl J Med. (2005) 352:987–96. doi: 10.1056/NEJMoa043330
50. Stupp R, Hegi ME, Mason WP, van den Bent MJ, Taphoorn MJ, Janzer RC, et al. Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial. Lancet Oncol. (2009) 10:459–66. doi: 10.1016/S1470-2045(09)70025-7
51. Singh N, Miner A, Hennis L, Mittal S. Mechanisms of temozolomide resistance in glioblastoma - a comprehensive review. Cancer Drug Resist. (2021) 4:17–43. doi: 10.20517/cdr.2020.79
52. Vredenburgh JJ, Desjardins A, Herndon JE 2nd, Dowell JM, Reardon DA, Quinn JA, et al. Phase II trial of bevacizumab and irinotecan in recurrent Malignant glioma. Clin Cancer Res. (2007) 13:1253–9. doi: 10.1158/1078-0432.CCR-06-2309
53. Gilbert MR, Dignam JJ, Armstrong TS, Wefel JS, Blumenthal DT, Vogelbaum MA, et al. A randomized trial of bevacizumab for newly diagnosed glioblastoma. N Engl J Med. (2014) 370:699–708. doi: 10.1056/NEJMoa1308573
54. Tamura R, Tanaka T, Miyake K, Tabei Y, Ohara K, Sampetrean O, et al. Histopathological investigation of glioblastomas resected under bevacizumab treatment. Oncotarget. (2016) 7:52423–35. doi: 10.18632/oncotarget.v7i32
55. Xie D, Pei Q, Li J, Wan X, Ye T. Emerging role of E2F family in cancer stem cells. Front Oncol. (2021) 11:723137. doi: 10.3389/fonc.2021.723137
56. Kim H, Lin Q, Yun Z. BRCA1 regulates the cancer stem cell fate of breast cancer cells in the context of hypoxia and histone deacetylase inhibitors. Sci Rep. (2019) 9:9702. doi: 10.1038/s41598-019-46210-y
57. Zhang M, Zhang Q, Bai J, Zhao Z, Zhang J. Transcriptome analysis revealed CENPF associated with glioma prognosis. Math Biosci Eng. (2021) 18:2077–96. doi: 10.3934/mbe.2021107
58. Yang H, Liu X, Zhu X, Zhang M, Wang Y, Ma M, et al. GINS1 promotes the proliferation and migration of glioma cells through USP15-mediated deubiquitination of TOP2A. iScience. (2022) 25:104952. doi: 10.1016/j.isci.2022.104952
59. Zhao Y, He J, Li Y, Lv S, Cui H. NUSAP1 potentiates chemoresistance in glioblastoma through its SAP domain to stabilize ATR. Signal Transduct Target Ther. (2020) 5:44. doi: 10.1038/s41392-020-0137-7
60. Genkai N, Homma J, Sano M, Tanaka R, Yamanaka R. Increased expression of pituitary tumor-transforming gene (PTTG)-1 is correlated with poor prognosis in glioma patients. Oncol Rep. (2006) 15:1569–74. doi: 10.3892/or
61. Guo L, Ding Z, Huang N, Huang Z, Zhang N, Xia Z. Forkhead Box M1 positively regulates UBE2C and protects glioma cells from autophagic death. Cell Cycle. (2017) 16:1705–18. doi: 10.1080/15384101.2017.1356507
62. Hu L, Li X, Liu Q, Xu J, Ge H, Wang Z, et al. UBE2S, a novel substrate of Akt1, associates with Ku70 and regulates DNA repair and glioblastoma multiforme resistance to chemotherapy. Oncogene. (2017) 36:1145–56. doi: 10.1038/onc.2016.281
63. Aquilanti E, Miller J, Santagata S, Cahill DP, Brastianos PK. Updates in prognostic markers for gliomas. Neuro Oncol. (2018) 20:vii17–26. doi: 10.1093/neuonc/noy158
64. Byron SA, Van Keuren-Jensen KR, Engelthaler DM, Carpten JD, Craig DW. Translating RNA sequencing into clinical diagnostics: opportunities and challenges. Nat Rev Genet. (2016) 17:257–71. doi: 10.1038/nrg.2016.10
65. Ye B, Li Z, Wang Q. A novel artificial intelligence network to assess the prognosis of gastrointestinal cancer to immunotherapy based on genetic mutation features. Front Immunol. (2024) 15:1428529. doi: 10.3389/fimmu.2024.1428529
66. Bian T, Zheng M, Jiang D, Liu J, Sun H, Li X, et al. Prognostic biomarker TUBA1C is correlated to immune cell infiltration in the tumor microenvironment of lung adenocarcinoma. Cancer Cell Int. (2021) 21:144. doi: 10.1186/s12935-021-01849-4
67. Albahde MAH, Zhang P, Zhang Q, Li G, Wang W. Upregulated expression of TUBA1C predicts poor prognosis and promotes oncogenesis in pancreatic ductal adenocarcinoma via regulating the cell cycle. Front Oncol. (2020) 10:49. doi: 10.3389/fonc.2020.00049
68. Wang J, Chen W, Wei W, Lou J. Oncogene TUBA1C promotes migration and proliferation in hepatocellular carcinoma and predicts a poor prognosis. Oncotarget. (2017) 8:96215–24. doi: 10.18632/oncotarget.v8i56
Keywords: cancer stem cells, spatial transcriptomics, single-cell RNA sequencing, prognostic signature, TUBA1C
Citation: Cao L, Lu X, Wang X, Wu H and Miao X (2024) From single-cell to spatial transcriptomics: decoding the glioma stem cell niche and its clinical implications. Front. Immunol. 15:1475235. doi: 10.3389/fimmu.2024.1475235
Received: 03 August 2024; Accepted: 28 August 2024;
Published: 17 September 2024.
Edited by:
Pengpeng Zhang, Nanjing Medical University, ChinaReviewed by:
Ge Zhang, The First Affiliated Hospital of Zhengzhou University, ChinaAimin Jiang, Fudan University, China
Jing Zhang, University of South Dakota, United States
Zhijia Xia, Ludwig Maximilian University of Munich, Germany
Copyright © 2024 Cao, Lu, Wang, Wu and Miao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hao Wu, MTU5OTYxNjA3MjBAMTYzLmNvbQ==; Xiaye Miao, bWlhb3hpYXllQDE2My5jb20=