Integration analysis of single-cell and spatial transcriptomics reveal the cellular heterogeneity landscape in glioblastoma and establish a polygenic risk model

Abstract

Background: Glioblastoma (GBM) is adults' most common and fatally malignant brain tumor. The heterogeneity is the leading cause of treatment failure. However, the relationship between cellular heterogeneity, tumor microenvironment, and GBM progression is still elusive.

Methods: Integrated analysis of single-cell RNA sequencing (scRNA-seq) and spatial transcriptome sequencing (stRNA-seq) of GBM were conducted to analyze the spatial tumor microenvironment. We investigated the subpopulation heterogeneity of malignant cells through gene set enrichment analyses, cell communications analyses, and pseudotime analyses. Significantly changed genes of the pseudotime analysis were screened to create a tumor progress-related gene risk score (TPRGRS) using Cox regression algorithms in the bulkRNA-sequencing(bulkRNA-seq) dataset. We combined the TPRGRS and clinical characteristics to predict the prognosis of patients with GBM. Furthermore, functional analysis was applied to uncover the underlying mechanisms of the TPRGRS.

Results: GBM cells were accurately charted to their spatial locations and uncovered their spatial colocalization. The malignant cells were divided into five clusters with transcriptional and functional heterogeneity, including unclassified malignant cells and astrocyte-like, mesenchymal-like, oligodendrocytes-progenitor-like, and neural-progenitor-like malignant cells. Cell-cell communications analysis in scRNA-seq and stRNA-seq identified ligand-receptor pairs of the CXCL, EGF, FGF, and MIF signaling pathways as bridges implying that tumor microenvironment may cause malignant cells' transcriptomic adaptability and disease progression. Pseudotime analysis showed the differentiation trajectory of GBM cells from proneural to mesenchymal transition and identified genes or pathways that affect cell differentiation. TPRGRS could successfully divide patients with GBM in three datasets into high- and low-risk groups, which was proved to be a prognostic factor independent of routine clinicopathological characteristics. Functional analysis revealed the TPRGRS associated with growth factor binding, cytokine activity, signaling receptor activator activity functions, and oncogenic pathways. Further analysis revealed the association of the TPRGRS with gene mutations and immunity in GBM. Finally, the external datasets and qRT-PCR verified high expressions of the TPRGRS mRNAs in GBM cells.

Conclusion: Our study provides novel insights into heterogeneity in GBM based on scRNA-seq and stRNA-seq data. Moreover, our study proposed a malignant cell transition-based TPRGRS through integrated analysis of bulkRNA-seq and scRNA-seq data, combined with the routine clinicopathological evaluation of tumors, which may provide more personalized drug regimens for GBM patients.

Keywords: ScRNA-seq; glioblastoma; heterogeneity; immune infiltration; spatial transcriptomics; tumor microenvironment; tumor progress-related gene risk score.

Figure 1

This study’s design and flowchart. *P < 0.05, ** P < 0.01, *** P < 0.001.

Figure 2

Single-cell RNA sequencing analysis of GBM. (A) The UMAP clustering map shows malignant cells, CD8 Tex, M1, Monocyte, and Oligodendrocytes representing the GBM cell types. (B) The proportion of malignant cells, CD8 Tex, M1, Monocyte, and Oligodendrocytes in GBM. (C) Differential gene expression analysis showing up- and downregulated genes across malignant cells, CD8 Tex, M1, Monocyte, and Oligodendrocytes. (Adjusted p-value: red for< 0.01, black for ≥ 0.01). (D) Bubble chart showing enrichment results of malignant cells based on GSEA analysis. (Normalized Enrichment Score: red for ≥ 0, blue for<0).

Figure 3

Intra-tissue heterogeneity of malignant cells in scRNA-seq. (A) The UMAP clustering map shows AC-like, MES-like, NPC-like, OPC-like, unclassified GBM malignant cells (B)The composition ratio of subtypes of malignant cells (AC-like, MES-like, NPC-like, OPC-like, and unclassified GBM malignant cells) in GBM patients. (C) Differential gene expression analysis showing up- and downregulated genes across subtypes of malignant cells (AC-like Malignant, MES-like Malignant, NPC-like Malignant, OPC-like Malignant, and unclassified GBM malignant cells). (Adjusted p-value: red for< 0.01, black for ≥ 0.01). (D) Various oncogenic hallmarks were enriched in AC-like Malignant, MES-like Malignant, NPC-like Malignant, OPC-like Malignant, and unclassified GBM malignant cells based on GSVA analysis.

Figure 4

Cell-cell interaction networks of scRNA-seq. (A) CD8 Tex communicates with AC-like Malignant, MES-like Malignant, M1, and Monocyte through CXCL signaling. (B) AC-like Malignant, NPC-like Malignant, OPC-like Malignant, and unclassified malignant cells communicate with M1 and Monocytes through EGF signaling. (C) AC-like Malignant, NPC-like Malignant cells communicate with Oligodendrocytes through FGF signaling. (D)CD8 Tex, AC-like Malignant, MES-like Malignant cells, M1, and Monocyte interact through MIF signaling.

Figure 5

Pseudotime analysis of malignant cells in GBM scRNA-seq. (A) Pseudotime is shown in single-cell trajectories of subtypes of GBM malignant cells that monocle2 inferred. Pseudotime is shown in a gradient from dark blue to light blue, indicating the onset of pseudotime. (B)The Pseudotime heatmap shows gene expression dynamics of significantly labeled genes. Genes (rows) are clustered into two modules, and cells (columns) are sorted according to pseudotime. (C) Ridge plots show changes in biological processes along the pseudotime. (D)Ridge plots show changes in hallmarks along the pseudotime.

Figure 6

Screening of prognosis-associated genes. (A) Correlations between error rate and classification trees. (B) The relative importance of prognosis-associated genes.

Figure 7

Validation of the TPRGRS model and performance analysis. (A) The overview of each patient’s survival status, risk rating distributions, and the expression of 15 TPRGRS genes in the TCGA dataset and two validation datasets (CGGA and REMBRANDT). (B) In three datasets, Kaplan-Meier survival curves revealed a shorter OS of the high-risk group than that of the low-risk group. (C) ROC curve analysis of the risk scores in three datasets, respectively.

Figure 8

Independent prognostic value of TPRGRS. Analysis of the Univariate and Multivariate Cox regression analysis revealed risk score was strongly associated with OS in the TCGA training dataset (A), the CGGA validation dataset (B), and the REMBRANDT validation dataset (C). Clinical traits analysis between the TPRGRS low-and high-risk groups (D). Pair-wise Fisher’s Exact test p-values (*P< 0.01, ·P< 0.05).

Figure 9

Functional analysis of the gene expression profile between the TPRGRS low- and high-risk groups. (A) GO enrichments analysis of upregulated genes in high-risk group in the TCGA dataset. (B) GSVA analysis showed enrichment scores of hallmark pathways in the high-risk group. (C) GSEA analysis reveals results for fifteen hallmark pathways associated with TPRGRS.

Figure 10

Landscape of somatic mutations in the TPRGRS low- and high-risk groups. (A, B) Waterfall plots displaying the somatic gene mutations of the high-(A) and low-(B) risk groups. Each column represents an individual patient. (C, D) The coincident and exclusive associations across mutated genes in high-(C) and low-(D) risk groups. The color or symbol of each cell represents the statistical significance of each pair of genes’ exclusivity or co-occurrence. Green represents mutual co-occurrence, and brown represents exclusive mutation. Pair-wise Fisher’s Exact test p-values (*P< 0.01, ·P< 0.05).

Figure 11

Immune landscape and drug prediction of high-risk groups. (A) The heatmap showing immune responses infiltration in the low- and high-risk groups. (B) ssGSEA showing higher enrichment scores of immunocytes-associated functions in the high-risk group. (C) A heatmap showing the differential expression of immune checkpoints in low- and high-risk groups. (D) The distribution of IC50 of four compounds of the TPRGRS low- and high-risk groups (*P< 0.05, ** P< 0.01, *** P < 0.001).