- Research article
- Open Access
Classification of glioma based on prognostic alternative splicing
BMC Medical Genomics volume 12, Article number: 165 (2019)
Previously developed classifications of glioma have provided enormous advantages for the diagnosis and treatment of glioma. Although the role of alternative splicing (AS) in cancer, especially in glioma, has been validated, a comprehensive analysis of AS in glioma has not yet been conducted. In this study, we aimed at classifying glioma based on prognostic AS.
Using the TCGA glioblastoma (GBM) and low-grade glioma (LGG) datasets, we analyzed prognostic splicing events. Consensus clustering analysis was conducted to classified glioma samples and correlation analysis was conducted to characterize regulatory network of splicing factors and splicing events.
We analyzed prognostic splicing events and proposed novel splicing classifications across pan-glioma samples (labeled pST1–7) and across GBM samples (labeled ST1–3). Distinct splicing profiles between GBM and LGG were observed, and the primary discriminator for the pan-glioma splicing classification was tumor grade. Subtype-specific splicing events were identified; one example is AS of zinc finger proteins, which is involved in glioma prognosis. Furthermore, correlation analysis of splicing factors and splicing events identified SNRPB and CELF2 as hub splicing factors that upregulated and downregulated oncogenic AS, respectively.
A comprehensive analysis of AS in glioma was conducted in this study, shedding new light on glioma heterogeneity and providing new insights into glioma diagnosis and treatment.
Glioma classification based on molecular characteristics plays an increasingly important role in diagnosis and treatment of glioma. Genetic, DNA methylation, gene expression features have been demonstrated to influence the prognosis of glioma patients [1, 2], and related molecular signatures, such as isocitrate dehydrogenase genes 1 and 2 (IDH1/IDH2) status, O6-methylguanine-DNA methyltransferase (MGMT) promotor status, codeletion of chromosome arm 1p and 19q (1p/19q codel) and TERT promoter status have been applied widely in prognosis prediction . Previous research on glioma classification has assisted in the prediction of prognosis and the guidance of treatment regimens, for example, the transcriptional classification of glioma reported by The Cancer Genome Atlas (TCGA) defined 4 subtypes, classical (CL), mesenchymal (MES), neural (NE), and proneural (PN), advancing our knowledge for the improvement of glioma diagnosis and therapy [1, 4,5,6]. However, glioma remains a serious threat to patients, especially the most aggressive kind, glioblastoma (GBM). Due to the lack of an effective treatment, the median survival of GBM patients is only 14.6 months after current standard therapy . We focused on a novel approach, classification of glioma based on alternative splicing (AS) event profiles, to understand glioma more comprehensively and explore new ideas for its diagnosis and treatment.
AS is a dynamic process that occurs after a gene is transcribed into precursor mRNA, and leads to complexity of the transcriptome. Changes in splicing patterns affect protein structure and function.
AS plays critical roles in oncogenesis. Cancer-specific mRNA transcripts, may result in loss-of-function of tumor suppressors or activation of oncogenes and cancer pathways . Most recently, a study comprehensively analyzed different tumor types from TCGA datasets to detect tumor-specific AS in combination with proteomic analysis; this study proposed a new method for exploring peptides as potential antigens in tumor immunotherapy . In addition to identifying tumor-specific AS, the biological effects of AS on cancer progression, especially prognosis-related AS, also warrant attention. Comprehensively analysis of non-small cell lung cancer , colorectal cancer  and esophageal carcinoma  revealed the prognostic significance of AS. AS of ANXA7, MARK4, MAX, USP5, WWOX, BIN, RON, and CCND1 were reported to affect critical biological functions of glioma, resulting in altered prognosis [13,14,15]; however, a systematic analysis of glioma splicing profiles has not been performed. Building on the availability of RNA-seq data of TCGA, we performed a comprehensive analysis of splicing events associated with prognosis in patients with glioma. Moreover, based on prognosis-related splicing events, we identified novel classification of glioma with distinct splicing characteristics and clinical features, shedding light on the new ideas for glioma research.
RNA-seq data and corresponding clinical data for 665 glioma samples (154 GBM, 511 LGG) were acquired from the data portal for TCGA (https://portal.gdc.cancer.gov/;DbGaP Study Accession:phs000178), using R/Bioconductor package TCGAbiolinks. Percent Spliced In (PSI) values, the percentage of splicing events in the abovementioned samples, were downloaded by using TCGASpliceSeq  (http://bioinformatics.mdanderson.org/TCGASpliceSeq), a web-based platform that provides splicing patterns of TCGA tumors. In this platform, splicing events were divided into 7 categories, including Exon Skip (ES), Retained Intron (RI), Alternate Promoter (AP), Alternate Terminator (AT), Alternate Donor site (AD), Alternate Acceptor site (AA), and Mutually Exclusive Exons (ME). In total, 83,776 splicing events were detected in GBM dataset, including 47,672 ES, 3574 RI, 12425 AP, 9210 AT, 4910 AD, 5630 AA and 355 ME. A total of 96,419 splicing events were detected in the LGG dataset, including 58,503 ES, 3694 RI, 13183 AP, 9366 AT, 5260 AD, 5954 AA and 458 ME.
The patients were divided into two groups by the mean cutoff PSI for each splicing event. To conduct further clustering analysis, which was needed to ensure that the data were not null and exclude the splicing events affected by outliers, splicing events that met the following conditions were included: 1) the PSI value was not missing for any samples; and 2) the sample size of each group was higher than 5% of the total size (n ≥ 8 in GBM, n ≥ 33 in glioma). Overall, 16,173 events across the pan-glioma samples and 20,939 events across the GBM samples matched. The log-rank test was used to assess the relationship between overall survival (OS) and AS events. Splicing events with a p-value< 0.01 in the log-rank test were regarded as prognostic. The same strategy was performed in additional studies, including survival analysis for splicing subgroups and splicing factors, as described in our former study [17,18,19].
Clustering analysis of prognosis-related splicing events was performed across all glioma (GBM + LGG) samples and for GBM samples.
Monte Carlo Consensus Clustering (M3C), a consensus clustering-based algorithm with a hypothesis testing framework, was used in our clustering study. The R package of M3C was downloaded from Bioconductor (https://www.bioconductor.org/). To evaluate the best number of subgroups (k), we comprehensively considered the following indexes: cumulative distribution function (CDF), empirical p value, proportion of ambiguous clustering (PAC) and relative cluster stability index (RCSI).
DNA methylation profiling
We evaluated the DNA methylation levels across splicing clusters by using a strategy similar to that reported by Michele Ceccarelli did .
Splicing events in each subtype were regarded as signatures if their mean PSI values were 30% higher than those in any other subtype (FDR < 0.05). Statistical analysis was performed by using the Mann-Whitney U test. The R package Scales was used to convert the PSI values of each splicing event to a scale from 0 to 1; subsequently, the converted values were displayed in heatmaps.
Correlation analysis of splicing factors and splicing events
Splicing factors with recognized prognostic significance in cancers, especially in glioma and pST1-specific splicing events, were selected for the analysis. A Spearman correlation test was used to evaluate the association between expression levels of splicing factors and PSI values of pST1-specific splicing events. Correlations with a coefficient > 0.45 and p-value< 0.05 were considered to be significant and are shown in the network diagram. Cytoscape (version 3.5.1) was used to construct the network of relationships.
IDH-status-related splicing events
Prognostic splicing events with mean PSI values that exceeded twice the difference between groups were regarded as IDH status-related splicing events (FDR < 0.05). Statistical analysis was performed by using the Mann-Whitney U test. The R package Scales was used to convert the PSI values of each splicing event to a scale from 0 to 1; subsequently, the converted values were displayed in heatmaps.
Identification of prognostic alternative splicing events in the TCGA GBM and LGG datasets
Splicing event profiles were analyzed in depth for 154 GBM patients and 511 LGG patients from TCGA (https://portal.gdc.cancer.gov/;DbGaP Study Accession:phs000178). Univariate survival tests were conducted to assess the correlation between the PSI of splicing events and OS (Additional file 1: Table S1, Additional file 2: Table S2). First, to search for prognostic splicing events in the pan-glioma samples, we combined the GBM and LGG datasets and conducted a comprehensive survival analysis; we identified 10,370 splicing events with prognostic significance, including 2611 ES, 693 RI, 1378 AP, 4456 AT, 568 AD, 632 AA and 32 ME (Fig. 1a, Additional file 1: Table S1). Further analysis revealed that there is a substantial difference in splicing patterns between GBM and LGG. Survival analysis was then performed for GBM alone. A total of 1038 splicing events were subsequently detected as significantly associated with GBM patient survival, including 344 ES, 45 RI, 164 AP, 322 AT, 85 AD, 77 AA and 1 ME (Fig. 1b, Additional file 2: Table S2). A total of 551 splicing events were involved in both pan-glioma and GBM prognosis (Fig. 1c). Interestingly, over 50% of splicing events were prognostic in each splicing category across pan-glioma samples, while less than 7% of splicing events were prognostic in each splicing category across GBM samples (Fig. 1d).
Consensus clustering identifies seven splicing types of glioma
To segregate the splicing subtypes across the pan-glioma samples, we performed M3C analysis for prognosis-related splicing events with 665 glioma samples. The data showed that the clustering stability was best at k = 2 or 7 (Fig. 2a, b, Additional file 7: Figure S1). Regardless of whether k = 2 or k = 7, glioma grade was the principal discriminator for GBM and LGG samples that were notably distributed in different clusters, demonstrating that clustering based on our criteria is consistent with the clinical classification of gliomas and indicating that splicing profiles differ between GBM and LGG (Fig. 2c). To explore the intrinsic characteristics of LGG, we further analyzed glioma splicing types at k = 7, labeled pST (pan-glioma Splicing Type) 1–7, which contained 153, 54, 173, 94, 20, 74, and 97 samples, respectively. Consistent with previous data, principal component analysis (PCA) based on signatures of pST1–7 in further analysis also revealed marked differences between GBM and LGG samples (Fig. 2d). As previously mentioned, pST1 was exclusively from the GBM dataset and included nearly all GBM samples (153/154, 99.4%) (Table 1). The other six clusters are almost all from LGG, except for one pST3 sample from GBM (Table 1).
To identify the clinical features and molecular characteristics of each cluster, we analyzed clinical data of each sample from TCGA and identified crucial molecular indicators of glioma. DNA methylation profiles were also shown to be clinically relevant for glioma classification . To determine the relationship between DNA methylation profiles and AS profiles, we performed an integrated analysis of data from TCGASpliceSeq and the TCGA methylation platform, including the HumanMethylation450/27 platforms.
Surprisingly, pST2–7 had distinct clinical and molecular features; among those types, the characteristics of pST2 were remarkable. Compared to pST3–7, pST2 contains a higher proportion of grade III gliomas (49/53, 92.5%, p < 0.01) and a higher proportion of astrocytoma (35/53, 66.0%, p < 0.01) (Table 1). Clinically, pST2 resembled pST1 or GBM, with a higher diagnosis age and an unfavorable prognosis (Table 1 and Fig. 2e). Interestingly, pST2 also shared characteristics similar to those of pST1 or GBM, with molecular enrichment for IDH wildtype (49/54, 90.7%, p < 0.01), 1p/19q non-codeletion (53/54, 98.1%, p < 0.01), MGMT promotor unmethylated (29/54, 53.7%, p < 0.01), chromosome 7 gain paired with chromosome 10 loss (37/54, 68.5%, p < 0.01), chromosome 19/20 gain (6/54, 11.1%, p < 0.01) (Table 1 and Fig. 3c). In addition, compared to pST3–7, pST1–2 showed genome-wide hypomethylation (Additional file 9: Figure S3). Overall, although the pST2 samples were all from the LGG dataset, this group had a malignant phenotype similar to that of GBM. Unlike pST2, pST7 harbored oligodendroglioma (58/84, 69.0%, p < 0.01), as indicated by histological analysis (Table 1 and Fig. 3c). Moreover, patients in pST7 had better clinical outcome, with molecular phenotypes characteristic of a favorable prognosis. Notably, pST7 was enriched for IDH mutation (96/97, 99.0%, p < 0.01), 1p/19q co-deletion (85/97, 87.6%, p < 0.01), and MGMT promotor methylation (94/97, 96.9%, p < 0.01), and it harbored no chromosome 7 gain paired with chromosome 10 loss or chromosome 19/20 gain (Table 1 and Fig. 3c).
Subsequently, we found our novel splicing clusters correlated to the published classifications of glioma, taking transcriptome subtype (CL, MES, NE, PN) for example, most of the pST1–2 samples were CL or MES (157/184, 85.3%), while pST6 were mainly NE (64/69, 92.8%) and pST3, 4, 5, 7 were mainly PN (p < 0.01) (Table 1 and Fig. 3c). Meanwhile, comparing to DNA methylation subtypes (LGm1–6), clear tendency of the distribution was observed: 92.6% (163/176) pST1–2 samples were divided into LGm4–6; while 90.2% (413/458) pST3–7 samples were divided into LGm1–3, which is genome-wide hyper-methylation (p < 0.01) (Table 1 and Fig. 3c).
Glioma classification analysis of genomics, RNA expression profiling, epigenetics or proteomics had been published. In this study, classification based on AS was also conducted. Actually, correlation among classifications were common , indicating that clinical phenomenon could be related to multiple molecular characteristics of different levels. Using these classifications comprehensively might be more reliable.
Consensus clustering identifies three splicing types of GBM
Compared to LGG, GBM has a distinct genotype and phenotype and a different splicing profile, according to our data. In addition, GBM requires more aggressive treatment due to its unfavorable outcomes and refractoriness. Heterogeneity in GBM is considered to contribute to treatment resistance [20, 21]. By using the single-cell RNA-seq technique, which was developed rapidly in recent years, intratumoral heterogeneity in GBM was identified at the single-cell level [22, 23]. To reveal the heterogeneity in GBM and offer more targeted suggestions for GBM diagnosis and treatment, classification of GBM according to molecular phenotype has always been the standard method [1, 2]. However, classification for GBM based on AS has not yet been explored. Using the same strategy used for pan-glioma classification, we performed a classification study for GBM alone. Monte Carlo Consensus Clustering of prognosis-related splicing profiles for 154 GBM samples identified three robust clusters (Fig. 3a, d), labeled ST (splicing type) 1–3, with the flattest curve in the CDF plot (Fig. 3b), the lowest p-value (Fig. 3c), and elbows in the PAC score curve and RCSI curve (Additional file 8: Figure S2A, S2B) at k = 3. ST1, ST2 and ST3 contain 49, 49 and 56 samples, respectively. PCA based on signatures of ST1–3 was also performed in further analysis to show the differences among the 3 splicing types of GBM samples (Fig. 3e). To identify the molecular characteristics of each cluster, we also analyzed several recognized molecular indicators for GBM, including IDH status, CpG island methylator phenotype (CIMP), and MGMT promoter status. In total, 6.8% (10/146) IDH mutations were observed in all GBM samples and differed among the three splicing clusters. IDH mutations occurred most frequently in ST3 (7/56, 12.5%), with no IDH mutations in ST1 (0/46) and a frequency of 6.3% (3/48) in ST2 (p < 0.05) (Table 2). Consistent with the association between IDH status and CIMP in glioma, ST1 was exclusively non-G-CIMP (48/48, 100%), while G-CIMP samples were enriched in ST2 (2/48, 4.2%) and ST3 (6/56, 10.7%) (p < 0.05) (Table 2). MGMT promoter status also differed among the three clusters. A higher proportion of MGMT promoter methylation samples was observed in ST3 (24/45, 53.3%) than in ST1 (11/37, 29.7%) or ST2 (17/41, 41.5%) (p = 0.099) (Table 2).
To further clarify the clinical significance of this novel classification based on survival-related splicing events, we analyzed the clinical data of TCGA patients. The most consistent clinical relation for previous GBM transcriptome classification was age, with a markedly higher proportion of younger patients in the proneural classification . In contrast, our GBM splicing subtypes were not significantly related to median age at diagnosis, but the proportion of younger patients (<=40) in ST1 (1/49, 2.0%) was smaller than that in non-ST1 (11/105, 10.5%) (Table 2). Compared to ST2 patients (15.1 months) and ST3 patients (15.6 months), ST1 patients (9.0 months) displayed markedly shorter median survival time (Table 2 and Fig. 3f). Although the difference was not statistically significant, longer disease-free survival was observed in ST3 patients (11.0 months) than in ST1 (8.4 months) and ST2 (8.5 months) patients (Table 2).
To assess the effect of standardized treatment on the 3 clusters, we examined TCGA data and compared the survival of patients receiving standardized treatment, defined as concurrent TMZ chemotherapy and radiotherapy, and nonstandard treatment. Standard treatment significantly improved prognosis in the ST3 cluster (p = 0.027), while it did not alter survival in the ST1 and ST2 clusters (Fig. 3g).
Moreover, we compared our classification to the previous transcriptional subtypes of GBM , and we found that ST1 samples were mainly distributed in MES (38/48, 79.2%) and rarely distributed in PN (1/48, 2.1%) (p < 0.01) (Table 2). Other comparisons of clustering by splicing profiles and clustering by previously reported classifications were also made to show the relationships between these classifications (Table 2).
In summary, clustering analysis of GBM based on prognostic AS events defined 3 subtypes: ST1 has the worst clinical outcome and molecular features of poor prognosis, ST3 has the best clinical outcome, treatment sensitivity and molecular features of favorable prognosis, and ST2 has clinical and molecular characteristics between those of the other 2 subtypes.
Identification of subtype-specific splicing events
To classify external glioma samples, subtype-specific splicing signatures were investigated. The splicing events that increased only in one subtype were identified and defined as signatures of each subtype. A heatmap was generated to display a total of 1175 splicing signatures (pST1–7: 265, 50, 5, 4, 370, 202, and 279) in each subtype of pan-glioma (Fig. 2c and Additional file 3: Table S3). The representative splicing signatures in each cluster were also presented in box plots, including AS of ZNF283 (AT), POLR2F (AT), LSM5 (ES), ZNF771 (AT), AKT3 (AT), GLS (AT), and LDHA (AP) (Fig. 4a). POLR2L, a subunit of RNA polymerase II, was found to be a prognostic splicing factor, and its AS was also prognostic in lung cancer . Here, we also discovered that another subunit of RNA polymerase II, POLR2L, whose AS was an important prognostic signature in pST2, might contribute to the malignant phenotype of low-grade glioma. AS of another splicing factor, LSM5 , was specifically upregulated in pST3. AS of AKT3 was reported to play important roles in cancers , and AKT3 expression was involved in glioma progression [26, 27]. The data demonstrated that AS of AKT3 was a splicing signature of pST5, which might further reveal the role of AKT3 in glioma. GLS and LDHA are metabolic enzymes that participate in biological progression of IDH-mutant glioma [28, 29]; however, the role of AS of GLS and LDHA in glioma remains unclear. Our data indicated that AS of GLS and LDHA were signatures of pST6 and pST7, respectively, the clusters with the best prognosis, suggesting that both expression and AS of metabolic enzymes might be involved in IDH-mutant glioma.
We performed the same strategy to identify 81 splicing signatures in GBM subtypes (ST1: 14, ST2: 4, ST3: 63) and draw a heatmap (Fig. 3d and Additional file 4: Table S4). Three representative splicing signatures were also shown in box plots (Fig. 4b). DHRS4L2 is a member of the dehydrogenase reductase family , with specific alternative transcripts expressed in neuroblastoma . AS of DHRS4L2 might also participate in GBM malignant progression because of its significant increase in ST1. In ST2, AS of FAM13A, a key regulator of lung cancer , was discovered as a splicing signature. RPL39L was reported as a ribosomal protein with upregulated expression in cancers and played a role in drug resistance [33,34,35]. Upregulated AS of RPL39L was observed in ST3 (Fig. 3g), moreover, PSI value of AS of RPL39L was negatively related to mRNA expression level of RPL39L (data not shown). These findings might indicate AS of RPL39L induce downregulation of mRNA expression, leading to sensitivity to concurrent treatment in GBM.
Regulatory network of splicing factors and pST1-specific splicing events
Changes in the expression levels of splicing factors significantly affect the development and progression of multiple tumors by regulating AS. Increasing evidence indicates that RNA processing triggered by alterations in splicing factors affects glioma outcomes . To elucidate prognostic AS regulatory processing in glioma, we characterized the regulatory network of pST1-specific splicing events, which represent unfavorable prognostic splicing signatures in glioma. Here, we focused on recognized prognostic splicing factors in tumors, especially in glioma, and identified their prognostic effect on glioma by conducting survival analysis using TCGA RNA-seq expression data. Splicing factors including cancer promoters (CLK2 , ELAVL1 , HNRNPA2B1 , HNRNPH1 , PTBP1 [13, 15, 40], SNRPB , SRSF1, SRSF2, SRSF3, SRSF7, SRSF9, SRSF10) and cancer suppressors (CELF2, MBLN2, QKI, RBFOX2 , RBM4, RBM5, RBM6, RBM10 , RBM11 , SRSF5) were applied to further correlation analysis. Significant correlations (p < 0.05, coefficient > 0.45 or < − 0.45) were shown in the network diagram (Fig. 5a and Additional file 5: Table S5), and 129 pST1-specific splicing events represented by their gene symbols (yellow nodes) are presented. Surprisingly, these pST1-specific splicing events were all positively (red lines) related to the expression of splicing factors associated with unfavorable prognosis (red nodes) and negatively (blue lines) related to the expression of splicing factors associated with favorable prognosis (blue nodes) (Fig. 5a). Meanwhile, we also found that expression level of above-mentioned splicing factors differed among splicing clusters, demonstrating that splicing factors play significant roles in splicing profiles of glioma (Additional file 10: Figure S4).
Recently, research revealed that SNRPB, the key element of spliceosome complex SmB/B′, is an important oncogenic splicing factor in GBM, contributing to the regulation of RNA processing, DNA repair, and chromatin remodeling . Consistent with this finding, the regulatory network of splicing factors and pST1-specific splicing events revealed SNRPB as a hub unfavorable prognostic splicing factor that is positively related to important pST1-specific splicing events (Fig. 5a). CELF2, a suppressor of glioma progression (Fig. 5b), was reported as an important regulator of many AS events across solid cancer samples  and exhibited a strong negative correlation with pST1-specific splicing events (Fig. 5a), indicating its key role in glioma suppression via the regulation of numerous splicing events. AS of KIF4A and AS of FKBP11 were hub events in this network (Fig. 5a) with important splicing predictors of glioma prognosis (Fig. 5c). The prominent correlations between expression of SNRPB and AT of KIF4A exon32 (coefficient = 0.70) and between expression of CELF2 and AT of FKBP11 exon8 (coefficient = − 0.72) are shown in scatter diagrams (Fig. 5d).
Identification of IDH-status-related splicing events
IDH status is the most significant indicator for prognosis prediction in glioma patients and was the predominant driver of transcriptome/methylome/fCNV glioma classification . In our splicing classification of pan-glioma, in addition to tumor grade, IDH status was also a vital discriminator of splicing clustering. The pST1–2 splicing subgroups harbored IDH-wild type samples (189/203, 93.1%), while pST3–7 harbored IDH mutations (411/455, 90.3%). To identify the IDH status-related prognostic splicing events, we separated the samples into IDH-wild type/mutant groups. Splicing events with mean PSI values that exceeded twice the difference between groups were chosen (FDR < 0.05), and a total of 840 splicing events were observed (Fig. 6a and Additional file 6: Table S6).
Although proteomics research indicated that a large number of genes have a predominant protein isoform , increasing evidences suggested that protein complexity derived from alternative splicing play a vital role in cancer progression [13, 48, 49]. AS may induce the insertion, deletion or substitution of functional domains in proteins, resulting in alteration of cellular function . To further predict the potential functional elements in coding proteins that are involved in glioma patient prognosis, we aligned the former data on IDH status-related splicing events with data on canonical protein products of genes in UniProt (www.uniprot.org), which contains annotations of known domains.
Zinc finger proteins (ZNFs), which contain ZNF domains, are an extensive family of proteins that contribute to varied biological functions, including cancer progression . Cys2His2(C2H2) was the first discovered and best characterized type of ZNF and has various functions, including sequence-specific DNA-binding and RNA binding . AS of ZNFs was observed to be related to IDH status. Moreover, the type of AS was basically AT, and the alternative exons were ZNF domains, namely, the alteration of deletion of the ZNF domains at the 3′ end correlated with IDH status. From the analysis results, a higher proportion of the deletions of the zinc finger motif in ZNFs (ZNF283, ZNF724P, ZSCAN20, ZNF606, ZNF169, ZNF430, ZNF20, KDM2B, etc.) were found in samples with IDH mutations (Fig. 6b, c), with similar results in pST2–7 samples (Fig. 4a), demonstrating that zinc fingers might contribute to malignant progression in glioma. Moreover, most of the alternative exons in these ZNFs encode tandem C2H2-type zinc finger motifs, except for KDM2B, which harbors CXXC-type and plant homeodomain (PHD)-type zinc finger motifs.
This study presents the first clustering analysis of prognosis-related AS in glioma. In contrast to previous glioma classification analysis of genomics, RNA expression profiling, epigenetics or proteomics [1, 2, 53], our study focused on AS, a biological process that affects the diversity of protein isoforms and functions.
Survival-related splicing events were identified systematically in pan-glioma and GBM samples from TCGA. A remarkably higher proportion of prognostic splicing events was observed across pan-glioma samples than in GBM samples. In addition, our data revealed clear differences between GBM and LGG samples. We conducted a clustering analysis of the merged GBM and LGG datasets and identified a total of 7 groups. The strong association between glioma grade (GBM versus LGG) and pan-glioma splicing classification that was observed in this study was in accordance with clustering based on proteomics . These results indicated that AS, as a co- or posttranscriptional process that affects protein translation, might contribute to glioma prognosis through influencing the histological phenotype and tumor grade.
As the only GBM-like cluster, pST1 had unique splicing characteristics compared to the remaining 6 LGG-like clusters, providing new methods to distinguish glioma grade. Furthermore, another significant finding of this study is that we identified pST2, an LGG cluster that harbored malignant phenotypes similar to those of GBM, predominantly IDH-wild type, indicating a novel method to distinguish low-grade gliomas with unfavorable prognoses based on iconic splicing events. The cluster-specific study also provided a new method to classify gliomas with different phenotypes using the detection of AS. Due to the limited samples in this study, larger cohorts are needed for verification.
Using the same strategy, the GBM cohort was separated into 3 STs. Among these types, ST1 was identified as the most distinctive subtype with the most unfavorable prognosis and was associated with recognized malignant phenotypic molecular characteristics. A previous GBM molecular classification study found that the proneural class is related to better outcomes but does not benefit from aggressive treatment . Conversely, our study suggested that ST3 had a favorable phenotype, with not only the best prognosis but also the strongest chemoradiotherapy sensitivity, indicating that splicing profiles may be independent of the transcriptome and affect GBM biological behavior. Furthermore, the study suggested that concurrent chemoradiotherapy might be particularly important for improving the prognosis of ST3 patients, but this possibility requires further clinical data.
Clustering of glioblastoma or pan-glioma based on splicing profiles had advantages in identifying clusters with clinically significant prognostic phenotypes, including malignant clusters, such as the ST1 GBM group and the pST2 LGG-like group, and benign phenotypes, such as the ST3 GBM group and the pST7 LGG-like group.
IDH status is considered to be the primary factor in most clustering studies , while a recent study revealed a link between IDH status and splicing events . Consistent with this finding, our data also demonstrated that prognostic splicing-based clustering was correlated with IDH status. However, the mechanism that underlies the relationship between IDH status and AS remains unclear. Although previous research revealed the correlation of spliceosome mutations and IDH status in primary myelofibrosis [54, 55], similar results were not observed in glioblastoma . One family that has a remarkable relationship with IDH status is ZNFs, which are also related to tumor grade; more deletions of ZNF domains were detected in IDH-mutant or LGG samples. Evidence showed that disruption of ZNF domains due to AS in ZNFs resulted in dysfunction of transcriptional regulation [57, 58]. Moreover, we found that most of the alteration of deletion of the ZNF domains correlated with IDH status was at the 3′ end, which might introduce premature stop codon (PTC) or nonsense mediated decay (NMD) into mRNA transcript, thereby regulate its stability . We further aligned the former AS analysis with corresponding RNA-seq data. Interestingly, negative relation between deletion of ZNF domain and mRNA expression level in some ZNFs (ZNF283, ZNF724P, ZSCAN20, ZNF606, ZNF20, etc.) were observed (Additional file 11: Figure S5), which might indicate this kind of AS introduced PTC or NMD and caused decay into mRNA transcript. However, this phenomenon was not universal in all ZNFs, demonstrating the underlying mechanism of deletion of ZNF domain in glioma was complicated and needed further investigation. Studies on the target genes of these alternative ZNF domains might provide insights into the significance of AS of ZNFs in glioma. Further studies that block the interaction between ZNFs and their target genes may reveal novel strategies to overcome glioma.
Regulatory network analysis of splicing factors and prognostic splicing events revealed a strong correlation between these two factors in glioma. PTBP1 and HNRNPs [13,14,15, 40] are recognized as important splicing factors in glioma progression that regulate various splicing events, which is also reflected in our analysis. We also found SNRPB and CELF2 as hub splicing factors of prognostic splicing events, indicating that these splicing factors might be potential targets for glioma treatment. Further research on the biological functions of these splicing factors and their downstream splicing events is needed.
In this study, a comprehensive analysis of AS in glioma was conducted. Novel classification of glioma based on AS profiles was established, which was associated with glioma grade. Regulatory network of splicing factors and prognostic splicing events were shown; suggesting hub splicing factors including SNRPB and CELF2 were participating in splicing regulatory in glioma. Moreover, we also identified IDH-status-related splicing events profiles. This study might shed new light on glioma heterogeneity and provide new insights into glioma diagnosis and treatment.
Availability of data and materials
RNA-seq data and corresponding clinical data for 665 glioma samples (154 GBM, 511 LGG) were acquired from the data portal for TCGA (https://portal.gdc.cancer.gov/;DbGaP Study Accession:phs000178), using R/Bioconductor package TCGAbiolinks, which are publicly available. Percent Spliced In (PSI) values, the percentage of splicing events in the abovementioned samples, were downloaded by using TCGASpliceSeq (http://bioinformatics.mdanderson.org/TCGASpliceSeq), a publicly available web-based platform that provides splicing patterns of TCGA tumors.
Alternate Acceptor site
Cumulative distribution function
CpG island methylator phenotype
Monte Carlo Consensus Clustering
Mutually Exclusive Exons
Proportion of ambiguous clustering
Percent Spliced In
Relative cluster stability index
The Cancer Genome Atlas
Verhaak RGW, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, Miller CR, Ding L, Golub T, Mesirov JP, et al. Integrated genomic analysis identifies clinically relevant subtypes of Glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell. 2010;17(1):98–110.
Ceccarelli M, Barthel FP, Malta TM, Sabedot TS, Salama SR, Murray BA, Morozova O, Newton Y, Radenbaugh A, Pagnotta SM, et al. Molecular profiling reveals biologically discrete subsets and pathways of progression in diffuse Glioma. CELL. 2016;164(3):550–63.
Eckel-Passow JE, Lachance DH, Molinaro AM, Walsh KM, Decker PA, Sicotte H, Pekmezci M, Rice T, Kosel ML, Smirnov IV, et al. Glioma groups based on 1p/19q, IDH, and TERT promoter mutations in tumors. N Engl J Med. 2015;372(26):2499–508.
Bhat K, Balasubramaniyan V, Vaillant B, Ezhilarasan R, Hummelink K, Hollingsworth F, Wani K, Heathcock L, James JD, Goodman LD, et al. Mesenchymal differentiation mediated by NF-kappaB promotes radiation resistance in glioblastoma. Cancer Cell. 2013;24(3):331–46.
Huse JT, Phillips HS, Brennan CW. Molecular subclassification of diffuse gliomas: seeing order in the chaos. GLIA. 2011;59(8):1190–9.
Ozawa T, Riester M, Cheng YK, Huse JT, Squatrito M, Helmy K, Charles N, Michor F, Holland EC. Most human non-GCIMP glioblastoma subtypes evolve from a common proneural-like precursor glioma. Cancer Cell. 2014;26(2):288–300.
Tanaka S, Louis DN, Curry WT, Batchelor TT, Dietrich J. Diagnostic and therapeutic avenues for glioblastoma: no longer a dead end? Nat Rev Clin Oncol. 2013;10(1):14–26.
Climente-González H, Porta-Pardo E, Godzik A, Eyras E. The functional impact of alternative splicing in Cancer. Cell Rep. 2017;20(9):2215–26.
Kahles A, Lehmann K, Toussaint NC, Hüser M, Stark SG, Sachsenberg T, Stegle O, Kohlbacher O, Sander C, Rätsch G, et al. Comprehensive analysis of alternative splicing across tumors from 8,705 patients. Cancer Cell. 2018;34(2):211–24.
Li Y, Sun N, Lu Z, Sun S, Huang J, Chen Z, He J. Prognostic alternative mRNA splicing signature in non-small cell lung cancer. Cancer Lett. 2017;393:40–51.
Xiong Y, Deng Y, Wang K, Zhou H, Zheng X, Si L, Fu Z. Profiles of alternative splicing in colorectal cancer and their clinical significance: a study based on large-scale sequencing data. EBIOMEDICINE. 2018;36:183–95.
Mao S, Li Y, Lu Z, Che Y, Sun S, Huang J, Lei Y, Wang X, Liu C. Zheng S, et al. CARCINOGENESIS: Survival-associated alternative splicing signatures in esophageal carcinoma; 2018.
Ferrarese R, Harsh GT, Yadav AK, Bug E, Maticzka D, Reichardt W, Dombrowski SM, Miller TE, Masilamani AP, Dai F, et al. Lineage-specific splicing of a brain-enriched alternative exon promotes glioblastoma progression. J Clin Invest. 2014;124(7):2861–76.
Golan-Gerstl R, Cohen M, Shilo A, Suh SS, Bakacs A, Coppola L, Karni R. Splicing factor hnRNP A2/B1 regulates tumor suppressor gene splicing and is an oncogenic driver in Glioblastoma. Cancer Res. 2011;71(13):4464–72.
Izaguirre DI, Zhu W, Hai T, Cheung HC, Krahe R, Cote GJ. PTBP1-dependent regulation of USP5 alternative RNA splicing plays a role in glioblastoma tumorigenesis. MOL CARCINOGEN. 2012;51(11):895–906.
Ryan M, Wong WC, Brown R, Akbani R, Su X, Broom B, Melott J, Weinstein J. TCGASpliceSeq a compendium of alternative mRNA splicing in cancer. Nucleic Acids Res. 2016;44(D1):D1018–22.
Li Y, Liu Y, Ren J, Deng S, Yi G, Guo M, Shu S, Zhao L, Peng Y, Qi S. miR-1268a regulates ABCC1 expression to mediate temozolomide resistance in glioblastoma. J Neuro-Oncol. 2018;138(3):499–508.
Yi G, Huang G, Guo M, Zhang XA, Wang H, Deng S, Li Y, Xiang W, Chen Z, Pan J, et al. Acquired temozolomide resistance in MGMT-deficient glioblastoma cells is associated with regulation of DNA repair by DHC2. BRAIN. 2019;142(8):2352–66.
Deng S, Li Y, Yi G, Lei B, Guo M, Xiang W, Chen Z, Liu Y, Qi S. Overexpression of COX7A2 is associated with a good prognosis in patients with glioma. J Neuro-Oncol. 2018;136(1):41–50.
Gilbert MR, Dignam JJ, Armstrong TS, Wefel JS, Blumenthal DT, Vogelbaum MA, Colman H, Chakravarti A, Pugh S, Won M, et al. A randomized trial of bevacizumab for newly diagnosed glioblastoma. N Engl J Med. 2014;370(8):699–708.
Nathanson DA, Gini B, Mottahedeh J, Visnyei K, Koga T, Gomez G, Eskin A, Hwang K, Wang J, Masui K, et al. Targeted therapy resistance mediated by dynamic regulation of extrachromosomal mutant EGFR DNA. SCIENCE. 2014;343(6166):72–6.
Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, Cahill DP, Nahed BV, Curry WT. Martuza RL, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. 2014;344(6190):1396–401.
Wang Q, Hu B, Hu X, Kim H, Squatrito M, Scarpace L, DeCarvalho AC, Lyu S, Li P, Li Y, et al. Tumor evolution of Glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. Cancer Cell. 2017;32(1):42–56.
Cui P, Zhang S, Ding F, Ali S, Xiong L. Dynamic regulation of genome-wide pre-mRNA splicing and stress tolerance by the Sm-like protein LSm5 in Arabidopsis. Genome Biol. 2014;15(1):R1.
Guo T, Sakai A, Afsari B, Considine M, Danilova L, Favorov AV, Yegnasubramanian S, Kelley DZ, Flam E, Ha PK, et al. A novel functional splice variant of AKT3 defined by analysis of alternative splice expression in HPV-positive Oropharyngeal cancers. Cancer Res. 2017;77(19):5248–58.
Joy A, Kapoor M, Georges J, Butler L, Chang Y, Li C, Crouch A, Smirnov I, Nakada M, Hepler J, et al. The role of AKT isoforms in glioblastoma: AKT3 delays tumor progression. J Neuro-Oncol. 2016;130(1):43–52.
Turner KM, Sun Y, Ji P, Granberg KJ, Bernard B, Hu L, Cogdell DE, Zhou X, Yli-Harja O, Nykter M, et al. Genomically amplified Akt3 activates DNA repair pathway and promotes glioma progression. Proc Natl Acad Sci. 2015;112(11):3421.
Chesnelong C, Chaumeil MM, Blough MD, Al-Najjar M, Stechishin OD, Chan JA, Pieper RO, Ronen SM, Weiss S, Luchman HA, et al. Lactate dehydrogenase A silencing in IDH mutant gliomas. Neuro-Oncology. 2014;16(5):686–95.
Panosyan EH, Lasky JL, Lin HJ, Lai A, Hai Y, Guo X, Quinn M, Nelson SF, Cloughesy TF, Nghiemphu PL. Clinical aggressiveness of malignant gliomas is linked to augmented metabolism of amino acids. J Neuro-Oncol. 2016;128(1):57–66.
Gabrielli F, Tofanelli S. Molecular and functional evolution of human DHRS2 and DHRS4 duplicated genes. GENE. 2012;511(2):461–9.
Zhang Q, Li Y, Liu G, Xu X, Song X, Liang B, Li R, Xie J, Du M, Xiao L, et al. Alternative transcription initiation and splicing variants of the DHRS4 gene cluster. Biosci Rep. 2009;29(1):47–56.
Eisenhut F, Heim L, Trump S, Mittler S, Sopel N, Andreev K, Ferrazzi F, Ekici AB, Rieker R, Springel R, et al. FAM13A is associated with non-small cell lung cancer (NSCLC) progression and controls tumor cell proliferation and survival. ONCOIMMUNOLOGY. 2017;6(1):e1256526.
Wong QW, Li J, Ng SR, Lim SG, Yang H, Vardy LA. RPL39L is an example of a recently evolved ribosomal protein paralog that shows highly specific tissue expression patterns and is upregulated in ESCs and HCC tumors. RNA Biol. 2014;11(1):33–41.
Ye Q, Ding SF, Wang ZA, Feng J, Tan WB. Influence of ribosomal protein L39-L in the drug resistance mechanisms of lacrimal gland adenoid cystic carcinoma cells. Asian Pac J Cancer Prev. 2014;15(12):4995–5000.
Dave B, Granados-Principal S, Zhu R, Benz S, Rabizadeh S, Soon-Shiong P, Yu K, Shao Z, Li X, Gilcrease M, et al. Targeting RPL39 and MLF2 reduces tumor initiation and metastasis in breast cancer by inhibiting nitric oxide synthase signaling. P NATL ACAD SCI USA. 2014;111(24):8838–43.
Marcelino Meliso F, Hubert CG, Favoretto Galante PA, Penalva LO. RNA processing as an alternative route to attack glioblastoma. Hum Genet. 2017;136(9):1129–41.
Park SY, Piao Y, Thomas C, Fuller GN, de Groot JF. Cdc2-like kinase 2 is a key regulator of the cell cycle via FOXO3a/p27 in glioblastoma. Oncotarget. 2016;7(18):26793.
Filippova N, Yang X, Wang Y, Gillespie GY, Langford C, King PH, Wheeler C, Nabors LB. The RNA-binding protein HuR promotes glioma growth and treatment resistance. Mol Cancer Res. 2011;9(5):648–59.
Lefave CV, Squatrito M, Vorlova S, Rocco GL, Brennan CW, Holland EC, Pan YX, Cartegni L. Splicing factor hnRNPH drives an oncogenic splicing switch in gliomas. EMBO J. 2011;30(19):4084–97.
Clower CV, Chatterjee D, Wang Z, Cantley LC, Vander HM, Krainer AR. The alternative splicing repressors hnRNP A1/A2 and PTB influence pyruvate kinase isoform expression and cell metabolism. Proc Natl Acad Sci U S A. 2010;107(5):1894–9.
Correa BR, de Araujo PR, Qiao M, Burns SC, Chen C, Schlegel R, Agarwal S, Galante PA, Penalva LO. Functional genomics analyses of RNA-binding proteins reveal the splicing regulator SNRPB as an oncogenic candidate in glioblastoma. Genome Biol. 2016;17(1):125.
Danan-Gotthold M, Golan-Gerstl R, Eisenberg E, Meir K, Karni R, Levanon EY. Identification of recurrent regulated alternative splicing events across human solid tumors. Nucleic Acids Res. 2015;43(10):5130–44.
Tsai YS, Dominguez D, Gomez SM, Wang Z. Transcriptome-wide identification and study of cancer-specific splicing events across multiple tumors. Oncotarget. 2015;6(9):6825–39.
Pavlyukov MS, Yu H, Bastola S, Minata M, Shender VO, Lee Y, Zhang S, Wang J, Komarova S, Wang J, et al. Apoptotic cell-derived extracellular vesicles promote malignancy of Glioblastoma via intercellular transfer of splicing factors. Cancer Cell. 2018;34(1):119–35.
Correa BR, de Araujo PR, Qiao M, Burns SC, Chen C, Schlegel R, Agarwal S, Galante PAF, Penalva LOF. Functional genomics analyses of RNA-binding proteins reveal the splicing regulator SNRPB as an oncogenic candidate in glioblastoma. GENOME BIOL. 2016;17(1).
Danan-Gotthold M, Golan-Gerstl R, Eisenberg E, Meir K, Karni R, Levanon EY. Identification of recurrent regulated alternative splicing events across human solid tumors. Nucleic Acids Res. 2015;43(10):5130–44.
Tress ML, Abascal F, Valencia A. Alternative splicing may not be the key to proteome complexity. Trends Biochem Sci. 2017;42(2):98–110.
Frampton GM, Ali SM, Rosenzweig M, Chmielecki J, Lu X, Bauer TM, Akimov M, Bufill JA, Lee C, Jentz D, et al. Activation of MET via diverse exon 14 splicing alterations occurs in multiple tumor types and confers clinical sensitivity to MET inhibitors. CANCER DISCOV. 2015;5(8):850–9.
Wan L, Yu W, Shen E, Sun W, Liu Y, Kong J, Wu Y, Han F, Zhang L, Yu T, et al. SRSF6-regulated alternative splicing that promotes tumour progression offers a therapy target for colorectal cancer. GUT. 2017:2017–314983.
Kornblihtt AR, Schor IE, Alló M, Dujardin G, Petrillo E, Muñoz MJ. Alternative splicing: a pivotal step between eukaryotic transcription and translation. NAT REV MOL CELL BIO. 2013;14(3):153–65.
Laity JH, Lee BM, Wright PE. Zinc finger proteins: new insights into structural and functional diversity. Curr Opin Struct Biol. 2001;11(1):39–46.
Pabo CO, Peisach E, Grant RA. Design and selection of novel Cys2His2 zinc finger proteins. Annu Rev Biochem. 2001;70:313–40.
Noushmehr H, Weisenberger DJ, Diefes K, Phillips HS, Pujara K, Berman BP, Pan F, Pelloski CE, Sulman EP, Bhat KP, et al. Identification of a CpG Island Methylator phenotype that defines a distinct subgroup of Glioma. Cancer Cell. 2010;17(5):510–22.
Lasho TL, Jimma T, Finke CM, Patnaik M, Hanson CA, Ketterling RP, Pardanani A, Tefferi A. SRSF2 mutations in primary myelofibrosis: significant clustering with IDH mutations and independent association with inferior overall and leukemia-free survival. BLOOD. 2012;120(20):4168–71.
Vannucchi AM, Lasho TL, Guglielmelli P, Biamonte F, Pardanani A, Pereira A, Finke C, Score J, Gangat N. Mannarelli C, et al. Mutations and prognosis in primary myelofibrosis. 2013;27(9):1861–9.
Brennan CW, Verhaak RGW, McKenna A, Campos B, Noushmehr H, Salama SR, Zheng S, Chakravarty D, Sanborn JZ, Berman SH, et al. The somatic genomic landscape of Glioblastoma. CELL. 2013;155(2):462–77.
Sayadi A, Jeyakani J, Seet SH, Wei CL, Bourque G, Bard FA, Jenkins NA, Copeland NG, Bard-Chapeau EA. Functional features of EVI1 and EVI1Delta324 isoforms of MECOM gene in genome-wide transcription regulation and oncogenicity. ONCOGENE. 2016;35(18):2311–21.
Liu Y, Huang W, Gao X, Kuang F. Regulation between two alternative splicing isoforms ZNF148(FL) and ZNF148(DeltaN), and their roles in the apoptosis and invasion of colorectal cancer. Pathol Res Pract. 2018.
Popp MW, Maquat LE. Nonsense-mediated mRNA decay and Cancer. CURR OPIN GENET DEV. 2018;48:44–50.
This study was supported by the National Natural Science Foundation of China (Grant nos. 81773290, 81472315, and 81872442), Guangdong Science and Technology Department (2017A030313497, 2016A040403053), PLA Logistics Research Project of China (CWH17L020, 18CXZ030). The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Prognostic Splicing Events in Pan-glioma.
Prognostic Splicing Events in GBM.
Specific Splicing Events of pST7.
Specific Splicing Events of ST3.
Correlation between Splicing Factors and pST1-specific Splicing Events (coeffcient > 0.45 or < − 0.45).
IDH-status-related splicing events.
Related to Fig. 2. A P-values for k=2 to k=8. B PAC score curve for k=2 to k=8. C RCSI curve for k=2 to k=8.
Related to Fig. 3. A PAC score curve for k=2 to k=8. B RCSI curve for k=2 to k=8.
DNA methylation among glioma splicing types. A DNA methylation among pST1-7. B DNA methylation among ST1-3.
Relative expression level of 22 splicing factors in pST1-7 glioma samples.
Scatter diagrams of PSI values of zinc finger proteins alternative splicing events and related gene expression levels.
About this article
Cite this article
Li, Y., Ren, Z., Peng, Y. et al. Classification of glioma based on prognostic alternative splicing. BMC Med Genomics 12, 165 (2019). https://doi.org/10.1186/s12920-019-0603-7
- Alternative splicing