Prognostic role of alternative splicing events in head and neck squamous cell carcinoma
Cancer Cell International volume 20, Article number: 168 (2020)
Aberrant alternative splicing (AS) is implicated in biological processes of cancer. This study aims to reveal prognostic AS events and signatures that may serve as prognostic predictors for head and neck squamous cell carcinoma (HNSCC).
Prognostic AS events in HNSCC were identified by univariate COX analysis. Prognostic signatures comprising prognostic AS events were constructed for prognosis prediction in patients with HNSCC. The correlation between the percent spliced in (PSI) values of AS events and the expression of splicing factors (SFs) was analyzed by Pearson correlation analysis. Gene functional annotation analysis was performed to reveal pathways in which prognostic AS is enriched.
A total of 27,611 AS events in 15,873 genes were observed, and there were 3433 AS events in 2624 genes significantly associated with overall survival (OS) for HNSCC. Moreover, we found that AS prognostic signatures could accurately predict HNSCC prognosis. SF-AS regulatory networks were constructed according to the correlation between PSI values of AS events and the expression levels of SFs.
Our study identified prognostic AS events and signatures. Furthermore, it established SF-AS networks in HNSCC that were valuable in predicting the prognosis of patients with HNSCC and elucidating the regulatory mechanisms underlying AS in HNSCC.
Recently, substantial progress in the area of high-throughput sequencing technology has motivated cancer genome research. Alternative splicing (AS) is a crucial posttranscriptional biological process that facilitates transcript variants and reprogramming of protein diversity in cells . It is known that dysregulation of AS often causes aberrant cellular homeostasis and is associated with malignancies [2, 3]. Accumulating evidence suggests that AS is implicated in the carcinogenesis and progression of cancers [4, 5]. Furthermore, researchers have found the clinical significance of AS, and it may serve as a prognostic predictor . For instance, Liu et al. reported that that DOCK5 variant promoted proliferation, migration, and invasion of HPV-negative HNSCC cells, and patients with higher expression of DOCK5 variant showed decreased overall survival . A recent study revealed that a novel splice variant of LOXL2 promotes progression of human papillomavirus-negative head and neck squamous cell carcinoma . These studies mainly focus on specific genes, however, studies providing a comprehensive evaluation of splicing events in HNSCC are scarce.
HNSCC is the sixth most common malignancy globally, and remains one of the leading causes of cancer-related death . HNSCC patients require a multimodal approach including surgical resection, radiotherapy, and systemic chemotherapy . However, about half of the patients will develop locoregional recurrence or metastasis, and their five-year survival rates remains unsatisfactory . These facts highlight the urgent need to identify underlying molecular mechanisms to develop effective therapy and improve the overall survival of patients.
RNA sequencing data generated by The Cancer Genome Atlas (TCGA) program enables researchers to illustrate the global profiling of AS events and identify the clinical significance and prognostic value of AS. TCGA SpliceSeq  provides valuable processed data for the analysis of AS events in 33 types of cancers, and it includes the following seven types of AS events: alternate acceptor site (AA), alternate donor site (AD), alternate promoter (AP), alternate terminator (AT), exon skip (ES), mutually exclusive exon (ME), and retained intron (RI).
In the present study, we attempted to elucidate the pattern of global aberrant AS events and its clinical and prognostic implications in patients with HNSCC using AS and clinical data obtained from TCGA database. Prognostic AS events that might function as prognostic indicators were identified and AS prognostic signatures were constructed. To assess the regulatory relationships between SFs and AS in HNSCC, a regulatory network was also established.
Materials and methods
TCGA data process
TCGA SpliceSeq is a database that provides AS profiles for seven types of splice events across 33 tumors based on TCGA RNA-seq data. To quantify AS events, the percent spliced in (PSI) value was processed for further data analysis. The PSI value indicates the inclusion of a transcript element divided by the total number of reads for that AS event. Alterations in PSI values range from 0 to 100%, which suggests a shift percentage in splicing events. AS events with a PSI value of more than 75% in a HNSCC cohort were obtained from TCGA SpliceSeq website. AS events with a standard diversion < 1 were excluded. An UpSet plot generated by the package “UpSetR”  in R software was used for analyzing and displaying the distribution and intersection among seven types of AS. Clinical information of patients with HNSCC was also downloaded and extracted from TCGA database. The overall survival (OS) was used as the endpoint for survival.
Survival analysis of AS events
In the survival analysis, the follow-up periods ranged from 91 to 6417 days after excluding patients with an OS of less than 90 days. Univariate Cox analysis was conducted to assess the correlation between the survival status of patients with HNSCC and PSI value (from 0 to 100) of each AS event (P < 0.05). A total of 486 patients with HNSCC were ranked from low to high according to risk score, and the median were used as cut-off to separate patients into low- or high- risk groups. Therefore, there are 243 patients in each group.
Prognostic signature construction
The top 20 most significant events of seven types of AS from the univariate Cox analysis were submitted to a least absolute shrinkage and selection operator (LASSO) analysis to develop prognostic signatures. The coefficients and partial likelihood deviance were calculated with the “glmnet”  package in R. The prognostic signatures for OS prediction were calculated by multiplying the PSI values of prognostic indicators with the coefficients assigned by LASSO Cox analysis. Evaluation of the splicing-based prognostic signature as an independent predictor was conducted by integrating the following clinical parameters into the multivariable Cox regression analysis: age, gender, tumor stage, lymph node status, distant metastasis, tumor-node-metastasis (TNM) stage and histological grade. The prognostic prediction efficacy of the AS signatures was examined by a time-dependent receiver-operator characteristic (ROC) curve and analyzed using the package “survivalROC”  in R software. The risk score calculated by the package “pheatmap”  was to evaluate the performance of the prognostic signatures. The Kaplan–Meier survival analysis was conducted to assess the survival difference between high- and low-risk groups.
SF-AS regulatory network
A list of 404 splicing factors was obtained from the database of SpliceAid 2 . The expression profiles of SF genes were selected from TCGA dataset. Correlation between the PSI values of prognostic AS events and the expression of SFs was examined by Pearson’s correlation analysis. SF-AS relationships with P value less than 0.05 and a Pearson correlation coefficient more than 0.65 were selected to construct the SF-AS a regulatory network via Cytoscape version 3.6.1.
Functional annotation of genes with prognostic AS events was performed by the package “clusterProfiler”  in R to investigate the functional relevance of the genes involved in AS events. Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Gene ontology (GO) were used to assess the functional categories. KEGG and GO terms with a P-value and q-value both smaller than 0.05 were considered significant categories.
Profiles of alternative splicing events in HNSCC
We processed TCGA splice-seq data and clinical information of a HNSCC cohort in TCGA, and a total of 347 patients were included in the analysis. In total, 27,611 AS events in 15,873 mRNAs were observed in HNSCC, indicating that AS events are common in the development of HNSCC. The specific numbers and percentages of events and corresponding genes in seven types of AS were showed in Fig. 1a. An UpSet graph was generated to analyze the intersection among seven types of AS and to display the distribution of spliced genes in different splicing types (Fig. 1b). We found that one gene may have multiple types of splicing events, and that ES was the most predominant type.
Prognostic AS events
To reveal the prognostic significance of AS events in patients with HNSCC, a univariate Cox analysis on AS events was performed. The specific numbers and percentages of prognosis-associated AS events and corresponding genes were showed in Fig. 1c. Additionally, one gene could present two or more AS events that were markedly related to the OS of patients with HNSCC. The UpSet plot demonstrated that ES was still the leading prognostic AS type, and that a gene could have up to three prognostic events (Fig. 1d).
Prognostic signatures for patients with HNSCC
The 20 most significant prognostic events of each of the seven AS types are shown in Fig. 2a–g. By using the LASSO Cox analysis, we developed seven types of prognostic signatures based on prognostic AA, AD, AP, AT, ES, ME, and RI events (Fig. 3a–g). Additionally, we conducted an integrated analysis of all the seven types of AS events to establish a comprehensive prognostic signature (abbreviated as “ALL”), which consist of RHOT1-40176-ES, SH3KBP1-88642-AP, AGTRAP-670-AA, SH3KBP1-88643-AP, PACS2-29633-AP, RBPMS-83289-AT, B3GNTL1-44424-AP, MOBP-64191-AT, NPHP3-66813-ES, ABCC5-67820-RI, FKTN-87134-ES and FKBP8-48446-AA (Fig. 3h). The detailed information of these eight prognostic signatures is listed in Table 1. As expected, Kaplan–Meier analysis suggested that the seven prognostic signatures could effectively separate the survival curves of high-risk groups from those of the low-risk groups (Fig. 4a–g), and the comprehensive prognostic signature could accurately predict prognosis (Fig. 4h). ROC curves further validated the efficacy of these eight prognostic signatures in prognosis prediction, and the area under the cure (AUC) of eight signatures was larger than 0.7, and the AUC of comprehensive signatures is 0.743 (Fig. 5a). We next performed univariate Cox regression analysis and found that the eight signatures had a high predictive value regarding the OS of patients with HNSCC (Fig. 5b). Furthermore, all eight signatures remained independent prognostic indicators for patients with HNSCC in multivariate analyses after other clinicopathological characteristics were adjusted (Fig. 5c–j). The risk scores of eight signatures in patients with HNSCC were ranked from low to high, as shown in Fig. 6a–h, and the median was used as a cut-off to divide high- and low-risk groups (upper panel). Patients with a high-risk score tended to have lower survival rate and shorter survival time (lower panel).
Prognostic SF-AS network
It is known that AS events are influenced by SFs . Hence, exploration of the SF-AS regulatory network is important to reveal the mechanism underlying AS in HNSCC. The Pearson correlation analysis indicated that there were 23 splicing factors positively correlated with 84 AS events, whereas 22 splicing factors were negatively correlated with 37 AS events. An interaction network was constructed according to the correlation between SF and AS, which comprises 28 SFs, 33 risk AS events (associated with poor prognosis), and 78 protective AS events (associated with good prognosis) (Fig. 7a). Among the network, splicing factors DDX39B, PRPF39, LUC7L3 and CLASRP were significantly correlated with more than 30 AS events. Notably, DDX39B directly regulates 86 AS events, therefore it was considered as a core SF.
Functional enrichment analysis
Functional enrichment analysis including the KEGG pathways enrichment and GO analysis were performed for survival‐related AS genes. The results of GO analysis indicated that AS genes were involved in biological processes such as autophagy and RNA splicing (Fig. 7b). In the KEGG analysis, genes corresponding to the prognosis-associated AS events were mainly enriched in cancer-related pathways, for instance, “autophagy”, “response to oxidative stress”, and “RNA splicing” (Fig. 7c).
AS is responsible for the modification of mRNA isoforms, and it plays an indispensable role in producing various mRNA and protein isoforms with multiple functions . Accumulating evidence has revealed that aberrant AS is implicated in the oncogenic processes of multiple malignancies [20, 21]. Therefore, investigation of AS mechanisms deepens our understanding of posttranscriptional regulatory patterns.
In recent years, next-generation sequencing technology has extensively promoted research at a whole-genome scale. RNA sequencing data from TCGA database have enabled the studies of AS patterns in various cancer types [22,23,24]. By using the SpliceSeq database, several studies explored alternative splicing profiles and constructed prognostic signatures for many types of cancers, including colorectal cancer , prostate adenocarcinoma , esophageal carcinoma , hepatocellular carcinoma [27, 28], kidney renal clear cell carcinoma , and soft tissue sarcoma . However, a comprehensive study regarding aberrant AS events in HNSCC is deficient. Li et al. carried out a systemic bioinformatic analysis on the genome-wide AS events of clinical HNSCC samples from the TCGA database . However, this study failed to identify AS events as independent prognostic predictors, and the authors did not explore the interaction between SFs and AS events. Another genome-wide analysis of the AS landscape in HNSC revealed novel AS events related to carcinogenesis and immune microenvironment, with implications for prognosis and therapeutic responses . This study revealed role of each individual AS events and genes in HNSCC immune microenvironment, instead of constructing of comprehensive AS prognostic signature or regulatory network. In the present study, a total of 27,611 AS events in 15,873 mRNAs were observed in HNSCC, indicating that AS events are common in the development of HNSCC. Results of the survival analysis suggest that 3433 AS events in 2624 genes are associated with the OS of patients with HNSCC. We constructed seven splicing prognostic signatures based on seven types of prognostic AS events. Additionally, a comprehensive prognostic signature was developed by integrating all seven types of AS. Genes and AS events enrolled in the comprehensive prognostic signature included RHOT1 (ES) , SH3KBP1 (AP) , AGTRAP (AA) , PACS2 (AP) , RBPMS (AT) , B3GNTL1 (AP) , MOBP (AT), NPHP3 (ES), ABCC5 (RI) , FKTN (ES) and FKBP8 (AA) , many of which are known to play important roles in cancer biology. For instance, RHOT1 can regulate cell migration and proliferation by suppressing the expression of SMAD4 in pancreatic cancer . PACS-2 as a phosphorylation-state dependent molecular switch that mediates either antiapoptotic or pro-apoptotic signaling . Mourskaia et al. suggest that ABCC5 functions as a mediator of breast cancer skeletal metastasis . FKBP8 is an endogenous inhibitor of mTOR, its degradation promotes tumor progression .
The AUC value of the comprehensive signature has reached 0.743, indicating that the prognostic biomarkers can be a useful tool to predict the prognosis of patients with HNSCC. These signatures could accurately distinguish between patients with HNSCC with distinct clinical outcomes, which confirmed that these prognostic signatures could serve as ideal predictors. Moreover, an SF-AS network was constructed to provide further insights into the regulatory mechanisms underlying AS in HNSCC. We found that DDX39B might act as core SFs because of their extensive correlation with AS events. Awasthi et al. found that increase in DDX39B enhances global translation and cell proliferation through upregulation of pre-ribosomal RNA, and dysregulation of DDX39B expression could lead to oncogenesis . Nakata identified DDX39 as a potential drug target for the treatment of AR splice variant-positive prostate cancer. The authors reported that DDX39B and its paralog DDX39A regulated androgen receptor splice variant AR-V7 generation .
HPV infection status is a critical factor in the carcinogenesis and progression of HNSCC. Therefore, we check the HPV status in the clinical data of TCGA-HNSC dataset as the reviewer suggested, and found that only 89 out of 528 patients were tested for HPV status (negative = 67, positive = 22). Subsequently, HPV positive and negative samples were integrated to identify differentially expressed alternative splicing. However, we failed to find significantly different AS events between HPV positive and negative samples with the threshold of |log2FC| > 1 and adj P Value < 0.05 (t-test). We think this is because the data of HPV status in most patients is absent, the differential analysis cannot provide accurate results.
We identified prognostic AS events based on the data from TCGA database. The constructed prognostic AS signatures could effectively predict survival outcomes of patients with HNSCC. The constructed prognostic SF-AS regulatory network may reveal the mechanisms underlying AS in the carcinogenesis of HNSCC. This in-depth analysis of AS in HNSCC may provide useful tools for prognosis prediction and clues for possible therapeutic targets for future clinical applications.
Availability of data and materials
Climente-González H, et al. The functional impact of alternative splicing in cancer. Cell Rep. 2017;20(9):2215–26.
Nilsen TW, et al. Expansion of the eukaryotic proteome by alternative splicing. Nature. 2010;463(7280):457–63.
El Marabti E, et al. The cancer spliceome: reprograming of alternative splicing in cancer. Front Mol Biosci. 2018;5:80.
Lee SC, et al. Therapeutic targeting of splicing in cancer. Nat Med. 2016;22(9):976–86.
Kozlovski I, et al. The role of RNA alternative splicing in regulating cancer metabolism. Hum Genet. 2017;136(9):1113–27.
Singh B, et al. The role of alternative splicing in cancer. Transcription. 2017;8(2):91–8.
Liu C, Guo T, Xu G, Sakai A, Ren S, Fukusumi T, Ando M, Sadat S, Saito Y, Khan Z, Fisch KM. Characterization of alternative splicing events in HPV-negative head and neck squamous cell carcinoma identifies an oncogenic DOCK5 variant. Clin Cancer Res. 2018;24(20):5123–32.
Liu C, et al. A novel splice variant of LOXL2 promotes progression of human papillomavirus-negative head and neck squamous cell carcinoma. Cancer. 2020;126(4):737–48.
Marur S, et al. Head and neck squamous cell carcinoma: update on epidemiology, diagnosis, and treatment. Mayo Clin Proc. 2016;91(3):386–96.
Budach V, et al. Novel prognostic clinical factors and biomarkers for outcome prediction in head and neck cancer: a systematic review. Lancet Oncol. 2019;20(6):e313–26.
Peitzsch C, et al. Cancer stem cells in head and neck squamous cell carcinoma: identification, characterization and clinical implications. Cancers. 2019;11(5):616.
Ryan M, et al. TCGASpliceSeq a compendium of alternative mRNA splicing in cancer. Nucleic Acids Res. 2016;44(D1):D1018–22.
Conway J, et al. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics (Oxford, England). 2017;33(18):2938–40.
Engebretsen S, et al. Statistical predictions with glmnet. Clin Epigenet. 2019;11(1):123.
Patrick JH, et al. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56(2):337–44.
Kolde R, et al. Robust rank aggregation for gene list integration and meta-analysis. Bioinformatics (Oxford, England). 2012;28(4):573–80.
Piva F, et al. SpliceAid 2: a database of human splicing factors expression data and RNA target motifs. Hum Mutat. 2012;33(1):81–5.
Sveen A, et al. Aberrant RNA splicing in cancer; expression changes and driver mutations of splicing factor genes. Oncogene. 2016;35(19):2413–27.
Yang X, et al. Widespread expansion of protein interaction capabilities by alternative splicing. Cell. 2016;164(4):805–17.
Wong ACH, et al. We skip to work: alternative splicing in normal and malignant myelopoiesis. Leukemia. 2018;32(5):1081–93.
Jin Y, et al. Mutually exclusive alternative splicing of pre-mRNAs. Wiley Interdiscip Rev RNA. 2018;9(3):e1468–e1468.
Xiong Y, et al. Profiles of alternative splicing in colorectal cancer and their clinical significance: a study based on large-scale sequencing data. EBioMedicine. 2018;36:183–95.
Lin P, et al. Systematic analysis of survival-associated alternative splicing signatures in gastrointestinal pan-adenocarcinomas. EBioMedicine. 2018;34:46–60.
Zhang Y, et al. Pan-cancer analysis of clinical relevance of alternative splicing events in 31 human cancers. Oncogene. 2019;38(40):6678–95.
Huang Z-G, et al. Prognostic value and potential function of splicing events in prostate adenocarcinoma. Int J Oncol. 2018;53(6):2473–87.
Mao S, et al. Survival-associated alternative splicing signatures in esophageal carcinoma. Carcinogenesis. 2019;40(1):121–30.
Chen Q-F, et al. Alternative splicing events are prognostic in hepatocellular carcinoma. Aging (Albany NY). 2019;11(13):4720–35.
Zhu G, et al. Prognostic alternative mRNA splicing signature in hepatocellular carcinoma: a study based on large-scale sequencing data. Carcinogenesis. 2019;40(9):1077–85.
Song J, et al. Systematic analysis of alternative splicing signature unveils prognostic predictor for kidney renal clear cell carcinoma. J Cell Physiol. 2019;234(12):22753–64.
Yang X, et al. Determining the prognostic significance of alternative splicing events in soft tissue sarcoma using data from The Cancer Genome Atlas. J Transl Med. 2019;17(1):283–283.
Li Z, et al. Systemic analysis of RNA alternative splicing signals related to the prognosis for head and neck squamous cell carcinoma. Front Oncol. 2020;10:87.
Li ZX, et al. Comprehensive characterization of the alternative splicing landscape in head and neck squamous cell carcinoma reveals novel events associated with tumorigenesis and the immune microenvironment. Theranostics. 2019;9(25):7648–65.
Li Q, et al. Role of RHOT1 on migration and proliferation of pancreatic cancer. Am J Cancer Res. 2015;5(4):1460–70.
Moss T, et al. Comprehensive genomic characterization of upper tract urothelial carcinoma. Eur Urol. 2017;72(4):641–9.
Zeng H, et al. Transcripto-based network analysis reveals a model of gene activation in tongue squamous cell carcinomas. Head Neck. 2019;41(12):4098–110.
Thomas G, et al. Caught in the act—protein adaptation and the expanding roles of the PACS proteins in tissue homeostasis and disease. J Cell Sci. 2017;130(11):1865–76.
Zhang Y, et al. Identification of a novel RBPMS-ROS1 fusion in an adolescent patient with microsatellite-instable advanced lung adenocarcinoma sensitive to crizotinib: a case report. Clin Lung Cancer. 2019;21(2):e78–83.
Sandanger T, et al. DNA methylation and associated gene expression in blood prior to lung cancer diagnosis in the Norwegian Women and Cancer cohort. Sci Rep. 2018;8(1):16714.
Mourskaia AA, et al. ABCC5 supports osteoclast formation and promotes breast cancer metastasis to bone. Breast Cancer Res. 2012;14(6):R149.
Hsu F, et al. Signal peptide peptidase promotes tumor progression via facilitating FKBP8 degradation. Oncogene. 2019;38(10):1688–701.
Awasthi S, et al. DDX39B promotes translation through regulation of pre-ribosomal RNA levels. RNA Biol. 2018;15(9):1157–66.
Nakata D, et al. The RNA helicase DDX39B and its paralog DDX39A regulate androgen receptor splice variant AR-V7 generation. Biochem Biophys Res Commun. 2017;483(1):271–6.
Ethics approval and consent to participate
The study was approved by the Ethics Committee of the Shaan Xi Provincial Tumor Hospital, and was performed following the TCGA publication guidelines.
Consent for publication
All authors have approved the submitted manuscript.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Ding, Y., Feng, G. & Yang, M. Prognostic role of alternative splicing events in head and neck squamous cell carcinoma. Cancer Cell Int 20, 168 (2020). https://doi.org/10.1186/s12935-020-01249-0