Glioblastoma cell differentiation trajectory predicts the immunotherapy response and overall survival of patients

Identification of 13 cell clusters in human GBMs using scRNA-seq data reveals high cell heterogeneity

A schematic diagram of the study design and principal findings is shown in Figure 1. Following the quality control standard and the normalization of GBM scRNA-seq data, 194 low-quality cells were excluded, and 2,149 cells from GBM cores were included in the analysis (Figure 2A). The number of genes detected was significantly related to the sequencing depth (Figure 2B). A total of 19,752 corresponding genes were included, and the variance analysis revealed 1,500 highly variable genes (Figure 2C). Principal component analysis (PCA) was performed to identify available dimensions and screen correlated genes. The top 20 significantly correlated genes are displayed as dot plots and heatmaps in Supplementary Figure 1. However, the PCA results did not demonstrate clear separations among cells in human GBMs (Figure 2D). We selected 20 principal components (PCs) with an estimated P value < 0.05 for subsequent analysis (Figure 2E).

Afterwards, the t-distributed stochastic neighbor embedding (tSNE) algorithm was applied, and cells in human GBMs were successfully classified into 13 separate clusters (Figure 2F). Differential expression analysis was performed, and a total of 8,025 marker genes from all 13 clusters were identified (Figure 2G). According to the expression patterns of the marker genes, these clusters were annotated by singleR and CellMarker (Figure 3A). Cluster 0, containing 518 cells, was annotated as GBM CSCs; clusters 1, 2, 6 and 10, containing 878 cells, were annotated as GBM cancer cells or GBM cells; cluster 3, containing 196 cells, was annotated as astrocytes; cluster 11, containing 44 cells, was annotated as oligodendrocytes; clusters 4, 5 and 9, containing 319 cells, were annotated as tumor-associated macrophages; cluster 8, containing 77 cells, was annotated as typical M1 macrophages; cluster 7, containing 81 cells, was annotated as typical M2 macrophages; and cluster 12, containing 36 cells, was annotated as T cells.

GBM cells can be divided into two subsets with distinct differentiation patterns, with their GDRGs correlating with immune and metabolic pathways

Trajectory analysis was performed to project all cells from GBMs onto one root and two branches, termed branches I and II (Figure 3B). The results showed that GBM CSCs were mainly located in the root, whereas GBM cells were mostly located in the branches. Interestingly, GBM cells in branch I, defined here as type I GBM cells, were all from clusters 2, 6, and 10 (434 cells), and GBM cancer cells in branch II, defined here as type II GBM cells, were all from cluster 1 (444 cells). The degree of differentiation of cells in the type I and type II GBM subsets varied significantly. The branch-dependent marker genes of type I and type II GBM cells were determined to be GDRGs. Differential analysis was performed, and 265 marker genes were identified as type I GDRGs, and 193 marker genes were identified as type II GDRGs. Therefore, a total of 498 GDRGs were ultimately identified in GBM (Supplementary Table 1).

Gene set enrichment analysis (GSEA) was performed to identify related molecular mechanisms and pathways of GBM cells with distinct differentiation patterns. The results revealed that type I GDRGs were significantly negatively correlated with the regulation of immune processes, such as antigen processing and immune cell differentiation (Figure 3C), and that type II GDRGs were significantly positively correlated with metabolism-related pathways, such as carbon metabolism, amino acid biosynthesis, glycolysis, and gluconeogenesis (Figure 3D). These findings from GDRGs indicate that GBM cells in distinct differentiation states demonstrate distinct tumor biology characteristics, which might provide new evidence for the molecular signatures of GBMs, including both intrinsic properties and the regulation of related pathways.

Then, we determined whether the observed GBM cell subsets could be identified using bulk RNA-seq data. As shown in Figure 4A–4C, the correlation analyses demonstrated that the type I GDRGs were highly correlated between scRNA-seq and bulk RNA-seq data from both the TCGA and CGGA databases, as were the type II GDRGs. These findings suggest that type I and type II GBM cells can also be identified by GDRG expression using bulk RNA-seq data, since these highly correlated genes could indicate a common cellular origin.

GDRGs in the two GBM cell subsets are functionally correlated and mostly mutated

To determine whether the gene profiles originating from different GBM cell subsets were functionally correlated, we utilized metagenes to represent the overall expression patterns of the corresponding gene profile. The expression of type I and type II metagenes, consisting of type I and type II GDRGs, respectively, were derived from the weighted averages of the expression of the constituent genes. As shown in Figure 4D–4F, the results demonstrated significantly strong correlations between type I and type II metagene expression in both scRNA-seq and bulk RNA-seq data, indicating that type I and type II GDRGs are functionally correlated even though they originate from distinct GBM subsets with different differentiation patterns.

We also analyzed the somatic mutation statuses of the GDRGs in the TCGA cohort. Most GDRGs (90.8%, 452/498) harbored mutations, and the top 9 mutated GDRGs with mutation frequencies ≥ 5% in GBM patients are shown in Figure 4G. Epidermal growth factor receptor (EGFR) exhibited the highest mutation frequency (53%), followed by CDK4 (16%) and TSPAN31 (16%). There were 246 genes (92.8%) exhibiting mutations among type I GDRGs and 262 genes (89.4%) exhibiting mutations among type II GDRGs. The mutation frequencies were not significantly different between the two groups (P=0.183, Figure 4H). However, 8 of the top 9 mutated GDRGs were type II GDRGs, and only 1 was a type I GDRG. These findings demonstrate the high mutation status heterogeneity in GDRGs, suggesting the pivotal roles of GDRGs in the development and progression of GBMs.

GDRG-based classification of GBM patients correlates with different OS outcomes and clinicopathological characteristics

To establish a classification of GBMs based on the expression patterns of GDRGs, machine learning-based unsupervised consensus clustering was performed on 151 GBM patients from the TCGA database. According to the relative change in the area under the cumulative distribution function (CDF) curve and consensus heatmap, the optimal number of clusters was determined to be two (k value = 2), and no appreciable increase was observed in the area under the CDF curve (Figure 5A–5C). Thus, all GBMs were classified into two groups: 80 (53.0%) in molecular cluster 1 (MC1) and 71 (47.0%) in molecular cluster 2 (MC2). Kaplan-Meier survival analysis indicated that patients with MC1 GBMs experienced significantly poorer OS than patients with MC2 GBMs (log-rank P=5.65×10^-3; Figure 5D).

Afterwards, we validated the GDRG-based classification in the CGGA cohort. As shown in Figure 5G–5I, the optimal number of clusters was also determined to be two (k value = 2), and the patients were also classified as MC1 (265 patients, 75.7%) and MC2 (85 patients, 24.3%). Kaplan-Meier survival analysis also demonstrated that patients with MC1 GBMs had poorer OS than patients with MC2 GBMs (log-rank P=1.26×10^-2; Figure 5J).

Additionally, we compared the expression patterns of the GDRGs and clinicopathological characteristics between two MCs of patients in the TCGA cohort. As shown in Figure 5E and 5F, the expression levels of the type I GDRG metagene were significantly higher (P<0.001) and the expression levels of the type II GDRG metagene were significantly lower (P<0.001) in patients with MC1 GBMs than in patients with MC2 GBMs. The same findings were also observed in the CGGA validation cohort (Figure 5K and 5L). Hence, we hypothesize that MC1 patients correspond mostly to the functional properties of the type I GBM subset and MC2 patients correspond mostly to the functional properties of the type II GBM subset.

The demographics and clinicopathological characteristics of patients with MC1 and MC2 GBMs are shown in Table 1. Compared with MC2 patients, MC1 patients were significantly younger in both the TCGA cohort (P=0.007) and CGGA cohort (P=0.002). However, no significant differences were observed in other variables between the two MCs of patients (all P > 0.05). Overall, the above findings indicate that this GDRG-based classification of GBM patients is robust and reliable across different populations, and different survival outcomes can be clearly discriminated according to this classification.

GDRG-based classification of GBM patients correlates with the expression patterns of immune checkpoints and different likelihoods of an immunotherapy response

The expression of 6 main immune checkpoints, namely, PDCD1 (PD1), CD274 (PDL1), PDCD1LG2 (PDL2), CTLA4, CD80 and CD86, was determined and compared between two GBM cell subsets and two MCs of GBM patients (Figure 6). In terms of scRNA-seq data, PD1, PDL1 and PDL2 were highly expressed in type I GBM cells, whereas CTLA4, CD80 and CD86 were highly expressed in type II GBM cells (Figure 6A). In terms of bulk RNA-seq data, PD1, PDL1 and PDL2 were highly expressed in MC1 GBM patients whereas CTLA4, CD80 and CD86 were highly expressed in MC2 GBM patients in both the TCGA database (Figure 6B) and CGGA database (Figure 6C).

Afterwards, the tumor immune dysfunction and exclusion (TIDE) algorithm was used to predict the likelihood of an immunotherapy response. According to the results from the TCGA training cohort, MC2 GBM patients (43.7%, 31/71) were more likely to respond to immunotherapy than MC1 patients (20.0%, 16/80, P = 0.003). Similarly, in the CGGA validation cohort, MC2 GBM patients (45.9%, 39/85) were more likely to respond to immunotherapy than MC1 patients (37.6%, 62/165, P < 0.001).

Then, SubMap analysis was performed to predict the likelihood of a clinical response to PD1 and CTLA4 inhibitors in the two MCs of GBM patients. MC1 GBM patients in both the TCGA and CGGA cohorts were more sensitive to anti-PD1 therapies, with Bonferroni-corrected P values of 0.001 and 0.033, respectively, while MC2 GBM patients were more sensitive to anti-CTLA4 therapies, with Bonferroni-corrected P values of 0.019 and 0.047, respectively (Figure 6D and 6E).

FN1, APOE, RPL7A and GSTM2 are the 4 most significant survival-predicting GDRGs in human GBMs

Univariate Cox analysis was performed, and 45 prognosis-associated GDRGs were identified in the TCGA training set. Least absolute shrinkage and selection operator (LASSO) followed by multivariate Cox analysis was then performed, and 4 significant survival-predicting GDRGs were identified (Supplementary Figure 2): fibronectin 1 (FN1, HR=5.28, P<0.001), apolipoprotein E (APOE, HR=0.39, P=0.001), ribosomal protein L7a (RPL7A, HR=0.27, P=0.021), and glutathione S-transferase mu 2 (GSTM2, HR=0.41, P=0.042). The expression levels of the 4 most significant prognostic GDRGs in both the scRNA-seq and bulk RNA-seq profiles are shown in Supplementary Figure 3. FN1 was significantly upregulated in GBM cells and T cells (general); APOE was significantly upregulated in GBM cells, CSCs and macrophages; RPL7A was significantly upregulated in GBM cells, CSCs and macrophages; and GSTM2 was significantly upregulated in GBM cells and astrocytes (Supplementary Figure 3A). Furthermore, the expression of these 4 GDRGs was validated using the Gene Expression Profiling Interactive Analysis (GEPIA) database, which includes 163 GBM (TCGA) samples and 207 normal (GTEx) samples. We found that all 4 survival-predicting GDRGs were upregulated in GBMs compared to normal tissues (Supplementary Figure 3B).

Kaplan-Meier survival analysis demonstrated that high expression of FN1 and low expression of APOE, GSTM2 and RPL7A were associated with poor OS in GBM patients (Supplementary Figure 3C).

Generation and validation of the GDRG-based prognostic risk score model to predict GBM patient survival

The prognostic risk score model was developed based on the above 4 GDRGs with the following formula: Risk score = Exp_FN1 × 1.66 + Exp_APOE × (-0.93) + Exp_RPL7A × (-1.30) + Exp_GSTM2 × (-0.90). Risk scores were calculated for all patients in the TCGA training set, and the patients were divided into either a high-risk (high score) group or a low-risk (low score) group using the median value of the risk score as the cutoff value (Figure 7C). Kaplan-Meier survival analysis demonstrated that patients in the high-risk group had significantly poorer OS than those in the low-risk group (log-rank, P = 2.778×10^-4; Figure 7A). The C-index of the GDRG signature for OS prediction was 0.781 (95% CI=0.742 to 0.820). Time-dependent receiver operating characteristic (ROC) analysis also indicated that the GDRG signature showed excellent performance in predicting the 0.5-, 1-, 2- and 3-year OS rates, with respective area under the curve (AUC) values of 0.767, 0.712, 0.752 and 0.776 (Figure 7B).

Then, the predictive formula was validated using the CGGA cohort in a similar manner. As shown in Figure 7F, all GBM patients were classified into high-risk or low-risk groups. Consistent with the results from the TCGA training set, the survival analysis also demonstrated that patients in the high-risk group had significantly poorer OS than patients in the low-risk group (log-rank, P = 9,783×10^-6; Figure 7D). The C-index of the GDRG signature was 0.715 (95% CI=0.676 to 0.754). Time-dependent ROC analysis also suggested favorable values in predicting OS in the CGGA validation set (Figure 7E). These results indicate that the GDRG-based prognostic risk score model can serve as a reliable prognostic predictor for different populations of GBM patients.

Development and validation of the clinically applicable prognostic nomogram with the GDRG signature and clinicopathological parameters

To investigate whether the prognostic significance of the GDRG signature is independent of the other clinicopathological variables in predicting the survival of GBM patients, univariate and multivariate Cox regression analyses were performed, and the results demonstrated that the GDRG signature was independently associated with OS in both the TCGA and CGGA cohorts (Table 2). Finally, a prognostic nomogram was successfully developed to provide a clinically applicable quantitative approach for individual OS prediction. Age, pharmacotherapy, radiotherapy, IDH mutation status, MGMT promoter methylation status, and the GDRG signature were included in the final OS prediction model (Figure 8A). The C-index of the prognostic nomogram was 0.896 (95% CI=0.857 to 0.935). Time-dependent ROC analysis revealed excellent predictive abilities for the 0.5-, 1-, 2- and 3-year OS rates, with AUC values of 0.734, 0.771, 0.864 and 0.919, respectively (Figure 8B). The calibration plots showed excellent agreement between the predicted 0.5-, 1- and 3-year OS rates and the actual observations in the TCGA cohort (Figure 8D–8F). Then, the prognostic model was validated in the CGGA cohort, with a C-index of 0.729 (95% CI=0.690 to 0.768). The time-dependent AUCs for the 0.5-, 1-, 2-, and 3-year OS rates with the prognostic nomogram were 0.725, 0.696, 0.694, and 0.701, respectively (Figure 8C). The calibration plots also showed excellent agreement between the OS predictions and the actual observations for the probabilities of the 0.5-, 1- and 3-year survival rates in the validation set (Figure 8G–8I). The above findings suggest the appreciable reliability of the prognostic nomogram for OS prediction that can be applied in different populations of GBM patients.

Data acquisition and preprocessing

The scRNA-seq data and bulk RNA-seq data of human GBM samples were included in this study for analysis. The scRNA-seq data of a total of 3,589 cells of 4 human primary GBM samples, accession number GSE84465 [34], were obtained from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database, containing 2,343 cells from tumor cores and 1,246 cells from peripheral regions, with a reading depth of 10× genomics based on Illumina NextSeq 500. The bulk RNA-seq profiles of GBM samples were obtained from the TCGA database (https://portal.gdc.cancer.gov/) and the CGGA database (http://www.cgga.org.cn). We excluded samples with unavailable clinical information and ultimately included 151 GBMs from the TCGA cohort as the training set and 350 GBMs from the CGGA cohort as the validation set.

Processing of the GBM scRNA-seq data

In total, 2,343 cells from tumor cores were included in this analysis. The Seurat package in R 3.5.1 was used for quality control, statistical analysis, and exploration of the scRNA-seq data [35]. First, 194 low-quality cells were excluded based on the following quality control standards: 1) genes detected in < 3 cells were excluded; 2) cells with < 50 total detected genes were excluded; and 3) cells with ≥ 5% of mitochondria-expressed genes were excluded. Then, the gene expression of the remaining 2,149 cells was normalized using a linear regression model. PCA was performed to identify significantly available dimensions with a P value < 0.05 [36]. Then, the t-distributed stochastic neighbor embedding (tSNE) algorithm was applied for dimensionality reduction with 20 initial PCs and for performing cluster classification analysis across all cells [37]. The differential expression analysis among all genes within cell clusters was performed using the limma package in R to identify the marker genes of each cluster. An adjusted P value < 0.05 and | log₂[fold change (FC)] | > 0.5 were considered the cutoff criteria for identifying marker genes. Afterwards, different cell clusters were determined and annotated by the singleR package according to the composition patterns of the marker genes and were then manually verified and corrected with the CellMarker database [38, 39]. The corresponding genes of cell surface markers for the annotation of cell clusters are listed in Supplementary Table 2.

Trajectory analysis and GDRG identification

Single-cell pseudotime trajectories of the GBM scRNA-seq data were constructed using the Monocle 2 algorithm [14]. This algorithm adopts a machine learning technique, learning a parsimonious principal graph to reduce the given high-dimensional expression profiles to a low-dimensional space. Single cells were projected onto this space and ordered into a trajectory with branch points. For data interpretation, the cells in the same branch were generally considered to be in the same differentiation state, while cells located in different branches were considered to have different cell differentiation characteristics. In addition, differential expression analysis was performed between branches, and genes that showed differential expression levels were defined as branch-dependent or state-specific genes or marker genes. These marker genes of GBM cells located in different branches were defined as GDRGs.

GSEA, correlation analysis, and somatic mutation analysis of branch-dependent GDRGs

GSEA (http://software.broadinstitute.org/gsea/index.jsp) was performed to identify the related molecular mechanisms and pathways of GBM cells in different differentiation states [40]. The adjusted P value was used to correct the false positive results by using the default Benjamini-Hochberg false discovery rate (FDR) method. An FDR value ≤ 0.05 of the enrichment gene sets was considered statistically significant. Then, we determined whether the observed cell subtypes generated from the scRNA-seq GBM data could also be identified in the bulk RNA-seq data. Because the coregulation of transcriptional programs in bulk RNA-seq data would indicate a common cellular origin, we believed that cell populations would be best identified by distinguishing gene expression profiles consisting of highly correlated genes [15]. Thus, correlation analyses of the expression of GDRGs were performed using Pearson’s correlation test, and the results were visualized as a heatmap. A P value < 0.05 and |correlation coefficient| > 0.3 were considered significantly correlative. To determine whether the gene profiles originating from different cell populations were functionally correlated, we utilized metagenes to represent the overall expression patterns of the corresponding gene profiles. The relative expression values of metagenes, an assembly of multiple genes, were derived as the weighted averages of the expression levels of the constituent genes using a logistic regression model [15]. Correlation analysis between branch-dependent GBM differentiation-related metagenes was performed by Pearson’s correlation test in both the scRNA-seq data and bulk RNA-seq data. Moreover, we performed somatic mutation analysis to identify the mutation statuses of the GDRGs. The cBioPortal database (https://www.cbioportal.org/) was used to download the somatic mutation data of GBM patients. The numbers of variant types and classifications were visualized with OncoPlot.

GDRG-based classifications of GBM patients in the TCGA and CGGA cohorts

Unsupervised consensus clustering, an algorithm based on k-means machine learning, was utilized to explore a molecular classification of both the TCGA and CGGA GBM cohorts based on the expression patterns of GDRGs using the ‘ConsensusClusterPlus’ package in R [41]. The clustering procedure, with 1,000 iterations, was performed by sampling 80% of the data in each iteration. The optimal number of clusters was determined by the relative change in the area under the CDF curves of the consensus score and consensus heatmap. Then, Kaplan-Meier survival analysis was performed to evaluate the prognosis of patients in different MCs. We also performed comparisons of the clinicopathological variables between different clusters of patients to further explore the associations between the GDRG-based MCs and the clinical features of GBM patients.

Immunotherapy response predictions

TIDE (http://tide.dfci.harvard.edu/) is a computational method that integrates the expression signatures of T cell dysfunction and exclusion to model tumor immune evasion [42]. We used the TIDE algorithm to predict the clinical response to immune checkpoint blockade (ICB) in GBM patients based on pretreatment genomics. The unsupervised subclass mapping method (SubMap; https://cloud.genepattern.org/gp/) was further applied to predict the response of GBM patients in different MCs to immunotherapy, including anti-PD1 and anti-CTLA4 therapy [43]. A Bonferroni-corrected P value < 0.05 was considered statistically significant.

Generation and validation of the GDRG-based prognostic risk score model

The associations between the expression levels of GDRGs and patient survival were first evaluated by univariate Cox regression analysis in the TCGA training cohort. The prognosis-related genes with a P value < 0.05 identified by the analysis were further screened by LASSO and multivariate Cox regression analyses. Consequently, a risk score model based on the key prognosis-related GDRGs was constructed to predict the prognosis of GBM patients [44]. We calculated the risk score of each patient by referring to our previously constructed formula, Risk score = Exp_GENE1 × β1 + Exp_GENE2 × β2 +…+ Exp_GENEn × βn, in which “Exp” represents the expression level of the corresponding gene, and “β” represents the regression coefficient calculated by the multivariate Cox analysis [9]. Accordingly, all GBM patients in the TCGA database were stratified into either the low-risk (low score) group or the high-risk (high score) group. Kaplan-Meier survival analysis was performed to estimate the OS of these two groups, and survival differences were evaluated by a two-sided log-rank test. The predictive accuracy of the GDRG-based prognostic model was evaluated by Harrell's concordance index (C-index) and time-dependent ROC curve analysis within 0.5, 1 and 3 years by utilizing the survcomp and survivalROC packages in R [33, 45]. Both the C-index and AUC ranged from 0.5 to 1, with 1 indicating perfect discrimination and 0.5 indicating no discrimination. Finally, the prognostic model generated by the TCGA training cohort was verified in the CGGA validation cohort.

Afterwards, univariate and multivariate Cox regression analyses were performed in both GBM cohorts to determine whether the predictive performance of the GDRG signatures could be independent of the clinicopathological variables. These variables included age, sex, Karnofsky performance status score, pharmacotherapy, radiotherapy, surgery, isocitrate dehydrogenase (IDH) mutation status, O6-methylguanine-DNA-methyltransferase (MGMT) promoter methylation status, telomerase reverse transcriptase (TERT) promoter mutation status, B-Raf proto-oncogene (BRAF) mutation status, X-linked alpha thalassemia mental retardation syndrome gene (ATRX) mutation status, EGFR mutation status, glioma CpG island methylator phenotype (G-CIMP) status, molecular subtype (classical, neural, mesenchymal, and proneural), and 1p/19q status.

Development and validation of the prognostic nomogram with the GDRG signature and clinicopathological parameters

Following univariate and multivariate Cox regression analyses, all the identified independent prognostic parameters were utilized to develop a prognostic nomogram for predicting the 0.5-, 1-, and 3-year survival outcomes of GBM patients using the rms package in R. Calibration plots at 0.5, 1, and 3 years were constructed to graphically evaluate the discriminative ability of the nomogram [45]. Then, the discrimination performance of the nomogram was quantitatively assessed by the C-index and ROC curve [33]. Finally, the nomogram constructed using the data from the TCGA cohort was validated in the CGGA cohort.

Data accessibility

The scRNA-seq data of GBM samples (accession number GSE84465) were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). The bulk RNA-seq profiles of GBM samples were obtained from the TCGA database (https://portal.gdc.cancer.gov/) and the CGGA database (http://www.cgga.org.cn).

Ethics approval

Ethics committee approval for our study was not required because the data were obtained from publicly available databases.