Integrated clustering signature of genomic heterogeneity, stemness and tumor microenvironment predicts glioma prognosis and immunotherapy response

Datasets and samples

The TCGA dataset (Supplementary Table 1) was downloaded from the Xena Browser at the University of California, Santa Cruz (UCSC, https://xenabrowser.net/datapages/) [28]. The clinical data (Supplementary Table 2) and RNA-seq data of the Chinese Glioma Genome Atlas (CGGA) were obtained from the CGGA data portal (http://www.cgga.org.cn/) [29]. The data for the IMvigor210 cohort were loaded from the R package “IMvigor210CoreBiologies” [30].

Analysis of genomic heterogeneity, stemness, and immune microenvironment

The glioma single nucleotide variant dataset processed by MuTect2 software [31] was downloaded from GDC (https://portal.gdc.cancer.gov/). Tumor mutation burden (TMB) was calculated for each glioma using the tmb function of the R package “maftools”, and microsatellite instability (MSI) for each glioma was obtained from previous studies [32, 33]. Glioma stemness scores based on RNA expression and DNA methylation were calculated according to a previous study [34]. StromalScore, ImmuneScore, and ESTIMATEScore were calculated for each glioma based on gene expression using the R software package “estimate” [35]. The abundance of tumor-infiltrating immune cells in gliomas was analyzed using the CIBERSOR, MCPCOUNTER and EPIC algorithms on the TIMER2 platform (http://timer.cistrome.org/) [36].

Consensus clustering analysis

Cluster analysis was performed by ConsensusClusterPlus [37], using agglomerative pam clustering with 1-pearson correlation distances, with 10 repetitions for 80% of the samples. The empirical cumulative distribution function plot was used to determine the optimal number of clusters.

Transfection with siRNA

These siRNAs were synthesized by the Shanghai GenePharma Co. (Supplementary Table 3). U87 MG cells were transfected with negative control siRNA and SH2D4A siRNA according to the manufacturer’s instructions of Lipofectamine 3000 (Invitrogen, L3000015). 48 hours after transfection, U87 MG cells were harvested for subsequent western blot and qRT-PCR experiments to verify their knock-down efficiency.

Western blot analysis

Western blot analysis was performed according to our previously reported method with minor modifications [38]. In brief, protein extracts were prepared using RIPA lysis buffer (Epizyme) in the presence of 100x EDTA-free Protease Inhibitor Cocktail (Epizyme). Cell lysates were separated through electrophoresis and electro-transferred to polyvinylidene difluoride (PVDF) membranes (Merck Millipore). The membrane was blocked with 5% non-fat milk in TBS solution for 2 h at room temperature and incubated with anti-SH2D4A (15957-1-AP; Proteintech) or anti-α-Tubulin (11224-1-AP; Proteintech) overnight at 4° C. Blots were then washed thrice with TBST and probed for 2 h at room temperature with HRP- conjugated AffiniPure Goat Anti-Rabbit IgG (SA00001-2; Proteintech). Finally, the labeled proteins were detected using the ECL reagent.

Quantitative real-time PCR

Trizol (Thermo Fisher) was used to extract RNA samples, and then 1 μg of total RNA and NovoScript® Plus All-in-one 1st Strand cDNA Synthesis SuperMix (gDNA Purge) (Novoprotein E047-01B) were used to prepare cDNA according to the manufacturer’s instructions. Thereafter, cDNA was then analyzed by NovoStart® SYBR qPCR SuperMix Plus (Novoprotein E096-01A). Relative gene expression was evaluated using the 2^−ΔΔCT method. The GAPDH was utilized as an internal control, and all primer sequences are compiled in Supplementary Table 4.

Cell proliferation and migration

Cell proliferation was analyzed by Cell Counting Kit-8 (CCK-8) (CellorLab, CX001L) and BeyoClick™ EdU-594 assay kit (Beyotime, C0078S). After 24h of transfection, cell suspensions were planted in 96-well plates. At 24, 48, 72 and 96h, 10ul CCK-8, and 90ul medium were incubated together at 37° C for 2h, then the absorbance at 450 nm was detected by a microplate reader.

For direct observation of the proliferating cells, the 5-Ethynyl-2-deoxyuridine incorporation experiment was also performed according to the specifications. The cells were further incubated with EdU for 2h before fixation, permeabilization, and EdU staining. The cell nuclei were stained with Hoechst33342 for 30min. Finally, the proportion of cells that incorporated EdU was detected using confocal laser microscopy.

To perform migration assays, U87 MG cells were seeded into 6-well plates. Transfection was performed when cell density reached 70-90%. 48 h later the plates were scraped with pipette tips, washed with PBS, and incubated with serum-free medium. The original images and migrated images were obtained using the inversion microscope system. The migrated area was analyzed by Image J software.

Construction of the risk model

In order to construct a scale-free co-expression network, a weighted correlation network analysis (WGCNA) was performed using the R package “WGCNA”, and four co-expression modules were finally obtained after merging modules with distances less than 0.25. Genes with high connectivity in the clinically important modules were identified as hub genes. LASSO and multivariate regression analyses were then performed sequentially to screen for positive hub genes that were significantly associated with overall survival (OS). Risk scores were calculated as follows:

$$Riskscore={\displaystyle \sum _{\text{i}=1}^{\text{n}}C{\text{oef}}_{\text{i}}\ast Ex{p}_i}$$

Identification of DEGs and enrichment analysis

We identified differentially expressed genes (DEGs) between different clusters by executing the “limma” package in the R software (p< 0.05 and |FC| ≥ 1.5). Gene enrichment analysis was conducted using Metascape [39].

Statistical analysis

One-way ANOVA, Wilcox test and t-test were used to analyze the significance of differences in heterogeneity, stemness, gene expression, and infiltration of immune cells in gliomas. LASSO, multivariate Cox regression, and Kaplan-Meier analyses were performed to screen and evaluate the risk signature using the R packages “glmnet” and “survival”. Roc curve was drawn using the R package “survivalROC.” All statistical analyses were performed using GraphPad Prism and R software, and P values less than 0.05 were considered statistically significant.

Availability of data and material

Bioinformatics datasets presented in this study can be found in online repositories, and the datasets used and/or analyzed during experiments are available from the corresponding author on reasonable request.

Consent for publication

All authors had final approval of the submitted versions and read the journal’s authorship statement.

Consensus clustering analysis of stemness, genomic heterogeneity, and the immune microenvironment in glioma

The molecular heterogeneity of gliomas significantly affects the overall prognosis of patients and their response to treatment. Three clusters were identified through unsupervised consensus clustering analysis of genomic heterogeneity of gliomas in the TCGA dataset (Figure 1A). We examined the overall survival of the three clusters and found significant differences in prognosis between them. Cluster 1 had the lowest prognosis, while cluster 3 had the highest prognosis (Figure 1B). Significant differences were observed in TMB and MSI between cluster 1, with the lowest prognosis, and cluster 3, with the highest prognosis (Figure 1C). In addition, gliomas were divided into two subtypes based on the unsupervised consensus clustering analysis of stemness scores in the TCGA dataset (Figure 1D). The prognosis of cluster 2 was more favorable than that of cluster 1 (Figure 1E). Assessment of glioma stemness scores based on RNA expression (RNAss) and DNA methylation (DNAss) revealed that cluster 1 had higher DNAss and lower RNAss than cluster 2 (Figure 1F). Tumor microenvironment is a crucial factor in glioma development and treatment resistance. We conducted another cluster analysis to examine the immune microenvironment of glioma. This analysis revealed that the StromalScore, ImmuneScore, and ESTIMATEScore differed significantly between clusters, with cluster 1 having higher scores than cluster 2. Improved prognosis was observed in cluster 2 (Figure 1G–1I).

Figure 1. Consensus clustering analysis in glioma. (A) Consensus clustering was performed based on the genomic heterogeneity of gliomas. (B) Kaplan-Meier curves displaying prognostic differences between different clusters. (C) The differences in TMB and MSI between clusters. (D) Consensus clustering was performed based on the stemness of gliomas. (E) Kaplan-Meier curves displaying prognostic differences between different clusters. (F) The differences in DNAss and RNAss between clusters. (G) Consensus clustering is performed based on the microenvironment of gliomas. (H) Kaplan-Meier curves displaying prognostic differences between different clusters. (I) The differences in the microenvironment scores between clusters.

Identification of DEGs associated with glioma stemness, genomic heterogeneity, and immune microenvironment

We quantified gene expression among the clusters divided according to genomic heterogeneity, stemness, and tumor microenvironment, and identified 1792 (412 down-regulated, 1380 up-regulated), 1599 (306 down-regulated, 1293 up-regulated), and 1806 (709 down-regulated, 1097 up-regulated) differentially expressed genes (DEGs), respectively (Figure 2A–2C). And finally, we identified 857 DEGs shared by the three clustering (Figure 2D). Enrichment analysis revealed that differentially expressed genes were involved in immunoregulatory pathways, including Cytokine-cytokine receptor interaction, Signaling by Interleukins, regulation of T cell activation, and Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell pathway (Figure 2E). We extracted all protein-protein interactions among differentially expressed genes from a PPI data source and constructed a PPI network. MCODE algorithm was applied to this network to identify neighborhoods where proteins are densely connected (Figure 2F). GO enrichment analysis was applied to the network to reveal biological meanings, and the results showed that MCODE 3 was enriched in HDACs deacetylate histones; MCODE 4 was enriched in TCR signaling Translocation of ZAP-70 to Immunological synapse, Phosphorylation of CD3 and TCR zeta chains; MCODE 6 was enriched in JAK-STAT signaling pathway, Signaling by Interleukins, and Cytokine Signaling in Immune system (Figure 2F).

Figure 2. DEGs screening and enrichment analysis. (A) Differentially expressed genes (Fold change >1.5 and P < 0.05) were screened between different clusters distinguished in genomic heterogeneity, (B) stemness, and (C) tumor microenvironment, respectively. (D) Venn diagram showing the overlap of identified differentially expressed genes. (E) Enrichment analysis of the differentially expressed genes. (F) Protein-Protein Interaction Networks, PPI.

Screening of risk hub genes

Through WGCNA analysis of integrated data containing DEGs and clinical characteristics (such as genomic heterogeneity, stemness and tumor immune microenvironment), four co-expression modules were finally determined (Figure 3A–3C). We filtered the most relevant module by evaluating its interrelationships with traits and identified the blue module as the most relevant (Figure 3D, 3E). The hub genes of blue module were extracted. Subsequently, LASSO Cox regression analysis and multivariate Cox regression analysis were performed sequentially to filter variables. Ultimately, we selected RAB42, SH2D4A, and GDF15 as risk hub genes (Figure 4A–4C). The expression of three risk hub genes in glioma was investigated in the TCGA and CGGA datasets, and RAB42, SH2D4A, and GDF15 were found to be higher expressed in GBM than in LGG (Figure 4D, 4E). In the TCGA and CGGA datasets, the Kaplan-Meier curve showed that RAB42, SH2D4A, and GDF15 significantly affected the prognosis of glioma, with a poor prognosis at high expression (Figure 4F, 4G).

Figure 3. Weighted correlation network analysis. (A, B) Soft-thresholding powers. (C) Clustering of module genes in the TCGA cohort. (D) Module-trait relationships. (E) Scatter plot of correlation between GS and MM. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Figure 4. Identification of risk hub genes related to the overall survival of glioma. (A, B) LASSO regression analysis of hub genes in the TCGA cohort. (C) Multivariate Cox regression analysis of hub genes in the TCGA cohort. (D, E) Violin plot showing the expression levels of the 3 risk hub genes between LGG and GBM in the TCGA and CGGA datasets, respectively. (F, G) Kaplan-Meier curves displayed 3 risk hub genes were significantly related to poor prognosis in the TCGA and CGGA datasets, respectively.

Construction of the risk signature

A risk model was then constructed in the TCGA dataset using the risk hub genes RAB42, SH2D4A, and GDF15. Kaplan-Meier curve showed that high-risk scores were associated with poor prognosis. The risk score had high sensitivity and specificity in predicting 1-, 3-, and 5-year survival rates in glioma patients (Figure 5A–5C). The results of the risk model were also confirmed in the CGGA dataset (Figure 5D–5F). In addition, we also validated the risk model in GBM and LGG, respectively (Supplementary Figure 1). The results of the correlation analysis showed that risk scores and risk hub genes RAB42, SH2D4A and GDF15 were highly correlated with the clinical traits such as TMB, MSI, RNAss, DNAss, StromalScore, ImmuneScore and ESTIMATEScore (Figure 6A). To better understand the relationship between the risk signature and clinical characteristics, we analyzed the distribution of survival status, WHO classification, risk scores, MGMT, and IDH status of glioma patients. We identified significant differences in clinical characteristics between patients with low-risk scores and those with high-risk scores. Patients in the low-risk score group exhibited an IDH-Mutant, Methylated phenotype and had a better prognosis, and GBM had a significantly higher risk score than the corresponding LGG subtype (Figure 6B–6E).

Figure 5. Construction of the risk score signature. (A) Distribution of the risk score, survival status, and expression profile of the prognostic genes in the TCGA cohort. (B) Kaplan-Meier curves displaying prognostic differences between high and low-risk groups in the TCGA cohort. (C) The ROC curves describing the sensitivity and specificity of the risk score in predicting OS at 1-, 3- and 5-year time points in the TCGA cohort. (D) Distribution of the risk score, survival status, and expression profile of the prognostic genes in the CGGA cohort. (E) Kaplan-Meier curves displaying prognostic differences between high and low-risk groups in the CGGA cohort. (F) The ROC curves describing the sensitivity and specificity of the risk score in predicting OS at 1-, 3- and 5-year time points in the CGGA cohort.

Figure 6. The relationship between risk score and clinical phenotype. (A) Correlations between 3 risk hub genes and tumor stemness, genomic heterogeneity, and immune microenvironment. (B) Sankey Diagram displayed the distribution of the survival status, IDH status, WHO grade, risk score and MGMT promoter status. (C) Analysis of the risk scores in different IDH status, (D) MGMT promoter status, and (E) WHO grades, respectively. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

The risk score is an independent risk factor for glioma prognosis

We performed univariate Cox regression analysis to investigate the independent prognostic factors for glioma. And the analysis showed that risk score, age, MGMT promoter status, WHO classification, and IDH status were significantly associated with prognosis (Figure 7A). Multivariate Cox regression analysis showed that risk score and age were significantly associated with prognosis, suggesting that risk score is an independent prognostic factor for glioma (Figure 7B). Following that, we developed a nomogram survival prediction model for glioma patients based on independent prognostic parameters, and plotted calibration curves. The results showed good agreement between the predicted outcome and the 1-, 3- and 5-year overall survival of patients (Figure 7C, 7D).

Figure 7. Risk score is an independent prognostic factor for glioma. (A) Univariate Cox regression analyses showed that clinical features such as the risk score, age, MGMT promoter status, WHO grade, and IDH status were significantly correlated with prognosis. (B) Multivariate Cox analysis showed the risk score remained associated with the prognosis. (C) The Nomogram was constructed to predict prognosis in patients at 1-, 3-, and 5 years in the TCGA dataset. (D) The calibration curve of the nomogram.

Risk signature is highly associated with tumor immune cells

The above analysis mentioned that DEGs were enriched in immunoregulatory process, and risk signature was highly correlated with the immune microenvironment, so we further analyzed the relationship between risk signature and infiltrating tumor immune cells. Infiltrating analysis showed elevated infiltration of T cell regulatory (Tregs), tumor associated macrophages, tumor associated fibroblasts, neutrophils and endothelial cells in patients in the high-risk score group compared to patients in the low-risk score group (Figure 8A–8C). Moreover, the risk signature and risk hub genes RAB42, SH2D4A, and GDF15 were highly positively correlated with the increase of these immunosuppressive tumor infiltrating cells (Figure 8D). Although CD4+ and CD8+ T cells increased in high-risk group, they might be exhausted T cells, characterized by progressive loss of T cell function and ultimately loss of cascade response. Therefore, we analyzed the relationship between risk signature and T cell exhaustion. The results showed that the risk signature was positively correlated with exhausted T cell signature (Figure 8E, 8F), and the expression of exhausted markers HAVCR2, TIGIT, LAG3, PDCD1, and LAYN was up-regulated in high-risk group (Figure 8G).

Figure 8. Risk signature is associated with infiltrating immune cells. (A) The level of immune cell infiltration between different risk subgroups evaluated by CIBERSORT, (B) MCPCOUNTER and (C) EPIC algorithm, respectively. (D) Correlation of infiltrating immune cells with risk scores and 3 risk hub genes. (E, F) Correlation of exhausted T cell signature with risk scores and 3 risk hub genes. (G) Violin plot showing the expression levels of the exhausted markers between low- and high-risk groups in the TCGA dataset. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Risk signature is associated with immunotherapy outcome

In order to analyze the predictive effect of risk signature on the efficacy of immune checkpoint inhibitors (ICI) therapy, we downloaded gene expression profiles and clinical data of the IMvigor210 cohort. All samples were divided into high and low-risk groups according to the risk signature. Patients with treatment response [complete remission (CR) or partial remission (PR)] had significantly lower risk scores than those without response [stable (SD) or progressive (PD)] (Figure 9A). By evaluating the distribution of CR/PR and SD/PD in the high-risk and low-risk groups, we found that the low-risk group responded better to ICI treatment than the high-risk group (Figure 9B), and had a significantly better prognosis than the high-risk group (Figure 9C). We then analyzed the prognostic value of risk signature and risk hub genes in pan-cancer (Figure 9D, 9E), and the results showed that the low-risk group had a better prognosis in LIHC\PAAD\MESO\UVM (Figure 9E).

Figure 9. Analysis of risk signature in immunotherapy cohort and pan cancer. (A) Violin plot depicting the risk scores between SD/PD and CR/PR groups in the IMVIgor210 cohort. (B) The response rate of the high-risk score group and low-risk score group to immune therapy in the IMVIgor210 cohort. (C) Survival analysis of patients with low-risk scores and high-risk scores in the IMVIgor210 cohort. (D) Survival analysis of 3 risk hub genes in pan cancer. (E) Kaplan-Meier survival curve of patients in high- and low-risk score groups in pan cancer. **, P < 0.01.

SH2D4A affects the migration and proliferation of glioma cells

To validate the expression of the risk hub genes in glioma, we analyzed the immunohistochemistry pathological specimen data. The results showed that the expression level of SH2D4A was increased in GBM relative to LGG (Figure 10A). However, the function of SH2D4A in glioma has not been reported in the literature so far. To investigate the effect of SH2D4A on glioma cells, we knocked down SH2D4A in U87 MG cells by transfecting specific siRNAs (Figure 10B, 10C). CCK-8 and Edu assays suggested that SH2D4A significantly affected the proliferation of U87 MG cells (Figure 10D–10F). Cell migration assays confirmed that knockdown of SH2D4A significantly inhibited the migration of U87 MG cells (Figure 10G–10H).

Figure 10. SH2D4A affects proliferation and migration of glioma cells. (A) Immunohistochemical staining of SH2D4A in glioma (data from HPA). (B, C) qRT-PCR and western blot analysis of SH2D4A knockdown efficiency in U87 MG cells. (D) Analysis of proliferation of control and SH2D4A-deficient U87 MG cells by CCK8 assay. (E, F) Representative images and statistical analysis of EdU assay in control and SH2D4A-deficient U87 MG cells. (G, H) Representative images and statistical analysis of cell migration assay in control and SH2D4A-deficient U87 MG cells at the indicated times. **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.