Jiting Qiu1, Chunhui Wang2, Hongkang Hu2, Sarah Chen3, Xuehua Ding2, Yu Cai1. 1. Department of Neurosurgery, Ruijin Hospital, Shanghai Jiaotong University School of Medicine, Shanghai 201803, P.R. China. 2. Department of Neurosurgery, Changzheng Hospital, Second Military Medical University, Shanghai 200003, P.R. China. 3. University of North Carolina at Chapel Hill School of Medicine, Chapel Hill, NC 27514, USA.
As the most aggressive primary central nervous system malignant tumour, glioblastoma (GBM) is a grade IV astrocytoma comprising 54% of all gliomas, with an incidence of 3.19 per 100,000 individuals in the United States in 2006–2010 (1). Although there are several therapeutic strategies for GBM treatment, including neurosurgical therapy, chemotherapy and radiotherapy, the median survival time for patients with GBM remains at 12–15 months, with a 5-year survival rate of <5% (1,2).At present, several molecular biomarkers, including epidermal growth factor receptor (EGFR), isocitrate dehydrogenase (IDH), O6-methylguanine DNA methyltransferase (MGMT) and PTEN have been tested in the clinical setting for patients with GBM (3). Given that the overall survival (OS) of patients with GBM remains low, novel molecular biomarkers and new treatment options are urgently required in order to determine the developmental mechanisms of GBM.With the rapid development of high-throughput sequencing and bioinformatics approaches, the study of oncogene expression has entered a new stage. To date, an increasing number of studies have proven that the results of genome-wide tumour bioinformatic analyses may be used as new biomarkers for diagnosis and treatment (4–7), which is important for investigating GBM-associated signaling pathways, such as the MAPK, PI3K and p53 signaling pathways. Multigene signatures have also been confirmed to predict the prognosis of patients with glioma based on mRNA expression profiling (8,9).Alternative splicing (AS) regulates the translation of mRNA isoforms and gives rise to protein diversity, thus serving as an important post-transcriptional regulatory mechanism (10). More than 95% of human genes undergo AS and encode splice variants in transcriptional processes (11). Increasing evidence indicates the essential role of AS in the course of oncogenesis, including tumour cell proliferation, immune escape, angiogenesis and tumour metastasis (12,13). In addition, specific splicing factor (SF) genes regulate AS events via binding to pre-mRNAs and yielding RNA splicing in the tumour microenvironment (14). Abnormal expression of SFs may result in the activation of oncogenes or the inactivation of cancer suppressors (15). Therefore, the role of specific SFs in pathogenesis provides theoretical support for tumour biological processes, especially at the gene transcription level (16). Constructing a prognostic model is essential to elucidate potential cancer biomarkers (17,18).Younger patients (median age, 36 years) with GBM generally have an improved prognosis and commonly carry IDH1 mutations, the cytosine-phosphate-guanine island methylator phenotype and a gene expression profile of a proneural subgroup. However, these biomarkers are valid for only a few primary patients with GBM (19). The present study constructed a prognostic model with good performance based on AS events for patients with GBM and plotted SF-AS networks that may serve as new molecular targets for the prognosis of patients with GBM.
Materials and methods
Data acquisition
Datasets used in the present study, including RNA sequencing (RNA-seq) data and corresponding clinical information of patients with GBM, were downloaded from The Cancer Genome Atlas (TCGA-GBM; http://tcga-data.nci.nih.gov/tcga/). All subtypes in this TCGA-GBM dataset, including classical, proneural, mesenchymal and neural, were analyzed without classification. The expression of SF genes in the mRNA splicing pathway was obtained from the SpliceAid2 database (http://www.introni.it/splicing.html). P<0.05 was considered to indicate a statistically significant difference.Different splicing types were classified using TCGA SpliceSeq (20), a Java application, to investigate the mRNA splicing patterns of RNA-seq and to identify significant changes in AS events. The percent-splice-in (PSI) value was calculated using the following formula with normalized read counts: , for seven common patterns of AS events, including alternative acceptor (AA), alternate donor (AD), alternate promoter (AP), alternate terminator (AT), exon skipping (ES), mutually exclusive exon (ME) and retained intron (RI) (21). In the current cohort, AS events with a PSI value >75% were obtained from the TCGA SpliceSeq database. The PSI value of AS events with standard deviation <1 were excluded from analysis.
Data analysis, dimension reduction and model construction
UpSet plot, a novel visualization technique for the quantitative results of multiple interactive sets (22), was used to visualize various combinations of the seven aforementioned AS types. To display the functional interactions of splicing-associated genes, a network was constructed using the Reactome FI plugin of Cytoscape (version 3.6.1) (23). This application predicts associations composed of specific genes and integrates them in a network plot.Univariate Cox regression was applied to analyze the association between AS events and OS to disclose the molecular characteristics of survival-associated AS events. Subsequently, the top 20 significant AS events of each type were used to develop prognostic predictor models.Least absolute shrinkage and selection operator (LASSO) Cox analysis, which is ideal for high-dimensional data (24), was performed to compute the optimal coefficient and the deviance likelihood for each prognostic feature using the ‘glmnet’ package in R (version 4.0–2; http://cran.r-project.org/web/packages/glmnet/index.html). According to each coefficient, the AS events were divided into high- and low-risk subgroups based on the median risk scores (All, 7.59; AA, 1.89; AD, 1.74; AP, 4.05; AT, 2.78; ES, 3.41; ME, 1.2; RI, 3.03). Kaplan-Meier survival analysis and log-rank test were used to further validate whether they resulted in diametrically distinct outcomes. Prognostic models were calculated by multiplying the PSI values of each significant splicing gene and the coefficient performed by LASSO Cox analysis.
Clinical prognostic analyses
Clinical parameters were obtained to assess the changes in AS events concerning the prognosis of patients with GBM. A total of 169 GBM samples and 5 normal tissues with available RNA-seq data were identified. Only cases with primary tumours, with no adjuvant pre-operative therapy and with ≥30 days of OS were included. There were 145 patients (51 females and 94 males; mean age, 59.83; age range, 21–85 years) with applicable clinical parameters and RNA-seq data who were included. To assess the efficiency of each prognostic candidate, the survivalROC package in R (version 1.0.3; http://cran.r-project.org/web/packages/survivalROC/index.html) was used to generate the area under the curve (AUC) of receiver operating characteristic (ROC) curves with censored data for each model (25).Evaluation of splicing-based prognostic signature as an independent predictor was performed by integrating the following clinical parameters into the univariate and multivariable Cox regression analysis: Age, sex, IDH status, MGMT promoter status and the risk score of AS events. The ability of the models to predict the survival outcome of patients with GBM was evaluated. All analyses were performed using R/Bioconductor (version 3.5.1; http://www.r-project.org/).
Correlation between splicing events and splicing factors
The correlation network between the gene expression levels of SFs and the PSI values of AS events were performed by Pearson's correlation analysis and plotted using Cytoscape (version 3.6.1; http://cytoscape.org/).
Results
Identification of survival-associated AS events
mRNA-seq datasets and clinical information of patients with GBM were obtained from TCGA (TCGA-GBM). A total of 169 GBM samples and 5 normal tissues with available RNA-seq data were identified. Only cases with primary tumours, with no adjuvant pre-operative therapy and with ≥30 days of OS were included. A total of 145 patients with applicable clinical parameters and RNA-seq data were finally included.For evaluation of prognostic values (Table I and Fig. 1A), there were a total of 3,827 alternate acceptor (AA) events in 2,684 genes, 3,269 AD events in 2,270 genes, 8,686 AP events in 3,476 genes, 8,456 AT events in 3,695 genes, 18,360 ES events in 6,935 genes, 184 ME events in 180 genes and 2,828 RI events in 1,897 genes. Hence, there may be ≥2 AS events in one gene associated with survival in patients with GBM.
Table I.
Summary of the glioblastoma sample cohort from The Cancer Genome Atlas.
(A) UpSet plots in glioblastoma, showing the RNA-sequencing distributions among the seven types of alternative splicing events. (B) Protein-protein interaction network with survival-associated splicing genes in glioblastoma. AA, alternate acceptor; AD, alternate donor; AP, alternate promoter; AT, alternate terminator; ES, exon skipping; ME, mutually exclusive exon; RI, retained intron.
Based on the univariate Cox regression analysis, a total of 115 AA events in 109 genes, 110 AD events in 106 genes, 346 AP events in 235 genes, 264 AT events in 179 genes, 631 ES events in 537 genes, 7 ME events in 7 genes and 96 RI events in 93 genes were identified as significant prognosis-associated AS events (P<0.05; Table I). Small nuclear ribonucleoprotein-associated protein N, heterogeneous nuclear ribonucleoprotein F, MRPS28, STAT3 and 26S proteasome non-ATPase regulatory subunit 4 were considered as hub genes in the network (Fig. 1B). Fig. 2A-G shows the top 20 significant survival-associated AS events based on PSI values. The two-sided red curves were obtained by the significant AS events in the volcano plot (Fig. 2H and Table SI).
Figure 2.
Top AS events associated with overall survival based on (A) AA, (B) AD, (C) AP, (D) AT, (E) ES, (F) ME and (G) RI. The larger and more red dots indicate the alternative splicing events with more significance. (H) Red dots represent splicing events that are significantly prognosis-associated (|z-score|>1). Blue dots represent splicing events without prognosis association (|z-score|<1). The x-axis of the z-score refers to either a positive or negative association. AS, alternative splicing; AA, alternate acceptor; AD, alternate donor; AP, alternate promoter; AT, alternate terminator; ES, exon skipping; ME, mutually exclusive exon; RI, retained intron.
Efficiency of prognostic models
Using LASSO Cox analysis, seven types of prognostic models were developed based on AA, AD, AP, AT, ES, ME and RI (Fig. 3). The risk score of each AS type was calculated with their PSI values, while high- and low-risk groups were divided using the median risk score as the cut-off point. The present study revealed that the survival time of the high-risk group was significantly shorter than that of the low-risk group in the current cohort (Fig. 4). Therefore, the prognostic models of each AS type were considered to predict the clinical outcome of patients with GBM (Fig. 4). In addition, the 16 most significant survival-associated AS events in the seven types were selected to construct the final prognostic model (Table II). The scatter plots and heat maps suggested that patients with high-risk scores had a low survival time, while patients with low-risk scores had a high survival time (Fig. S1). The final prognostic model was deemed to be an ideal predictor of what could significantly distinguish patients with GBM with distinct survival times (Fig. S1). The AUC of the ROC curve validated the performance of prognostic models with good performance in prognosis prediction (Fig. 5). Additionally, univariate Cox regression analysis was performed to assess the prognostic value of clinical parameters, including age, sex, IDH mutation, MGMT promoter methylation and risk score of AS events. Multivariate Cox regression was applied after the sex parameter was eliminated since there was no significant association between sex and survival in the univariate analysis. The hazard ratios for risk score of AS events were 1.0071 (95% CI, 1.0032–1.0111) and 1.0063 (95% CI, 1.0024–1.0103) by univariate and multivariate analyses, respectively (Table III). In the present cohort, the parameters of age, MGMT promoter methylation status and risk score of AS events were considered as independent factors of prognosis prediction.
Figure 3.
Construction of prognostic signatures based on least absolute shrinkage and selection operator Cox analysis. Each colored curve refers to a significant splicing event. While the log(λ) values of the lower horizontal coordinate increase, the coefficients of each splicing event tend toward stationarity. The values of the upper horizontal coordinate represent non-zero coefficients. AA, alternate acceptor; AD, alternate donor; AP, alternate promoter; AT, alternate terminator; ES, exon skipping; ME, mutually exclusive exon; RI, retained intron.
Figure 4.
Kaplan-Meier curves of prognostic models for high- and low-risk subgroups of patients with glioblastoma. AA, alternate acceptor; AD, alternate donor; AP, alternate promoter; AT, alternate terminator; ES, exon skipping; ME, mutually exclusive exon; RI, retained intron.
Table II.
Prognostic signatures based on each type of alternative splicing event.
AA, alternate acceptor; AD, alternate donor; AP, alternate promoter; AT, alternate terminator; ES, exon skip; ME, mutually exclusive exons; RI, retained intron; AUC, area under the curve.
Figure 5.
AUC of ROC curves of eight types of prognostic models for glioblastoma within 1, 2 and 3 years. AUC, area under the curve; ROC, receiver operating characteristic; AA, alternate acceptor; AD, alternate donor; AP, alternate promoter; AT, alternate terminator; ES, exon skipping; ME, mutually exclusive exon; RI, retained intron.
Table III.
Cox regression analysis of clinical parameters and risk score for assessing prognostic model value in patients with glioblastoma.
Univariate Cox regression
Multivariate Cox regression
Clinical variable
HR (95% CI)
P-value
HR (95% CI)
P-value
Age, >60 vs. ≤60 years
1.0397 (1.0210–1.0587)
2.57×10−05
1.0344 (1.0125–1.0567)
0.00189
Sex, male vs. female
1.0198 (0.6541–1.5899)
0.93117
N/A
N/A
IDH mutation, yes vs. no
0.1836 (0.0573–0.5881)
0.00432
0.4370 (0.1272–1.5019)
0.18876
MGMT promoter methylation, yes vs. no
0.4987 (0.3129–0.7950)
0.00345
0.5253 (0.3276–0.8425)
0.00755
Risk score of AS events, high vs. low
1.0071 (1.0032–1.0111)
0.00036
1.0063 (1.0024–1.0103)
0.00153
HR, hazard ratio; IDH, isocitrate dehydrogenase; MGMT, O6-methylguanine DNA methyltransferase; N/A, not applicable; AS, alternative splicing.
SF-AS regulatory association
The regulatory network of the significant survival-associated AS events (n=241) and SFs was plotted using Cytoscape (n=25) (Table SII and Fig. 6A). Among all SFs, four representative prognostic factors were selected, including heat shock protein b-1 (HSPB1), protein arginine N-methyltransferase 5 (PRMT5), protein FAM50B (FAM50B) and endoplasmic reticulum chaperone BiP (HSPA5). The expression levels of each SF gene were calculated, while high- and low-risk groups were divided using the median values as cut-off. The results revealed that the survival time of the high-risk group was significantly shorter than that of the low-risk group (P<0.05; Fig. 6B-E). Therefore, HSPB1, PRMT5, FAM50B and HSPA5 were identified as significantly representative prognostic factors. Among the splicing correlation network, a total of 129 favorable prognosis AS events were correlated with survival-associated SFs, while a total of 112 adverse prognosis AS events were correlated with survival-associated SFs (P<1×10−10; Table SIII). Notably, the most favorable splicing events were negatively regulated by SFs, while the most adverse splicing events were positively regulated by SFs (Table SIII).
Figure 6.
(A) Alternative splicing events whose percent-splice-in values are positively or negatively associated with overall survival are represented with green or red dots, respectively. Survival-associated splicing factor genes are represented with purple dots. The positive or negative correlation between splicing factor genes and splicing events are represented with red or green lines, respectively. Kaplan-Meier survival analysis was performed for (B) HSPB1, (C) PRMT5, (D) FAM50B and (E) HSPA5 as significantly representative prognostic factors. HSPB1, heat shock protein b-1; PRMT5, protein arginine N-methyltransferase 5; FAM50B, protein FAM50B; HSPA5, endoplasmic reticulum chaperone BiP.
Discussion
GBM is a primary neuroepithelial tumour of the central nervous system and accounts for 12–15% of all intracranial tumours (1,2). The present study analyzed GBM datasets composed of classical, proneural, mesenchymal and neural subtypes. Patients with GBM have a poor median survival time of 12–15 months following standard therapy, with only 3–5% of patients surviving up to 5 years after the first diagnosis (1,2). Currently, several molecular markers have been tested as part of the routine clinical investigation of patients with GBM, including MGMT, IDH, EGFR, PTEN, VEGF, TP53, p16INK4a gene and 1p19q gene, as well as imaging biomarkers (26). However, there is still a limited number of molecular signatures for the contribution to anti-GBM therapies, such as temozolomide, bevacizumab and lomustine (27). Developments in next-generation sequencing methods have led to the identification of specific molecular signatures of GBM that allow for further investigation of the molecular pathogenesis of this disease (28). In recent years, high-throughput RNA-seq approaches have extensively promoted genome-wide analyses, including genome splicing investigation. The present study used bioinformatics techniques to identify survival-associated AS events in order to construct splicing signatures for the prediction of prognosis, orchestrate SF-AS networks and assess their potential underlying molecular mechanism.Previously, SpliceSeq analyses have been adopted to establish AS profiling and construct prognostic models in glioma; several potential AS events were identified in pan-glioma and GBM cohorts, including adenine DNA glycosylase, metalloreductase STEAP3, SUMO-conjugating enzyme UBC9, von Hippel-Lindau disease tumor suppressor, BTB/POZ domain-containing protein KCTD7, protein S100-A4, endothelin-converting enzyme 2 and lymphocyte antigen 6K (29–31). Additionally, several prognostic models based on AS events have been constructed for prognosis prediction, which may complement the molecular classification, further identify potential glioma subgroups and highlight SFs as an important mechanism of splicing regulation in the carcinogenesis and aggressiveness of GBM (9,31–33).The present study comprehensively analyzed the prognostic value of AS events and SFs in a GBM cohort using several computational approaches. The in-depth study further investigated alterations of mRNA-seq for prognostic monitoring. The ideal prognostic model built by combining all significant AS events exhibited potential for predicting the survival outcome of patients with GBM. Splicing correlation network analysis further revealed regulated nodes, revealing the potential mechanisms in the regulatory network at the genome-wide level.In the interaction network analysis, HSPB1, PRMT5, FAM50B and HSPA5 were identified as independent prognostic factors. It has been reported that HSPB1 (also known as Hsp27) phosphorylation leads to the activation of orphan nuclear receptor TAK1 and TAK1-p38/ERK pro-survival signaling, thus acting against TNF-α-induced apoptosis (34). PRMT5 is one of the candidate genes required for apoptosis or loss of self-renewal for differentiated and undifferentiated GBM cells, respectively (35). The specificity and efficacy of four novel PRMT5 inhibitors have been identified for the treatment of GBM (36). Additionally, a previous study has validated that a family with sequence similarity to FAM50B serves a key role as methylation-based biomarkers for the diagnosis and treatment of GBM (37). Furthermore, by specifically inhibiting HSPA5, a new compound known as HA15 was able to increase the unfolded protein response and lead to the death of cancer cells by concomitant induction of autophagy and apoptosis, both in vitro and in vivo (38). Whether downregulation of specific SFs may affect the associated AS events requires further validation in vivo.Several limitations inevitably influenced the reliability of the present results. There was a limited number of patients with GBMs with complete clinicopathological parameters recruited in the present analysis. All subtypes, including classical, proneural, mesenchymal and neural, were analyzed without precise classification. Therefore, subsequent functional experiments in vitro and in vivo are required to further validate the molecular mechanisms of how SFs regulate the splicing process in GBM development.Finally, the present study identified that survival-associated AS events were favorable predictors and the prognostic model performed well in predicting the stratification for patients with GBM. According to these identified survival-associated AS events and SFs, several valuable biomarkers may be determined for further validation studies.In conclusion, the present study established a molecular phenomenon of OS-associated AS and SFs in patients with GBM, which is valuable for investigating the underlying mechanisms in the oncogenesis of GBM. The present findings may facilitate the ongoing effort in developing novel transcriptome prognostic models for the management of GBM. Further identification of prognostic SFs and construction of an SF-AS network will advance the investigation of splicing-associated mechanisms.
Authors: Y K Banasavadi-Siddegowda; L Russell; E Frair; V A Karkhanis; T Relation; J Y Yoo; J Zhang; S Sif; J Imitola; R Baiocchi; B Kaur Journal: Oncogene Date: 2016-06-13 Impact factor: 9.867
Authors: Alexander Lex; Nils Gehlenborg; Hendrik Strobelt; Romain Vuillemot; Hanspeter Pfister Journal: IEEE Trans Vis Comput Graph Date: 2014-12 Impact factor: 4.579
Authors: Wei Zhao; Jun Li; Mei-Ju M Chen; Yikai Luo; Zhenlin Ju; Nicole K Nesser; Katie Johnson-Camacho; Christopher T Boniface; Yancey Lawrence; Nupur T Pande; Michael A Davies; Meenhard Herlyn; Taru Muranen; Ioannis K Zervantonakis; Erika von Euw; Andre Schultz; Shwetha V Kumar; Anil Korkut; Paul T Spellman; Rehan Akbani; Dennis J Slamon; Joe W Gray; Joan S Brugge; Yiling Lu; Gordon B Mills; Han Liang Journal: Cancer Cell Date: 2020-11-05 Impact factor: 31.743
Figure 1.
Figure 2.
Figure 3.
Figure 4.
Figure 5.
Figure 6.