Systematic analysis and prediction model construction of alternative splicing events in hepatocellular carcinoma: a study on the basis of large-scale spliceseq data from The Cancer Genome Atlas

TCGA-based data collection

The genes expression profile data as well as the clinical information of HCC patients had been collected from TCGA (https://cancergenome.nih.gov/). Information of seven types AS events were collected from the TCGA SpliceSeq database (https://bioinformatics.mdanderson.org/TCGASpliceSeq/), including alternate acceptor site (AA), alternate promoter (AP), alternate donor site (AD), alternate terminator (AT), exon skip (ES), mutually exclusive exons (ME), and retained intron (RI). Due to defining perioperative mortality as death that occurs within 30 days of surgery may underestimate ‘true’ mortality among patients undergoing hepatic resection (Mayo et al., 2011). We only included patients with a survival time more than 90 days, finally, 335 HCC patients were selected in this study. The Percent Spliced In (PSI) value, rating from 0 to 1 which is used to quantify AS events in general. To generate a reliable set of AS events, we applied a series of stringent filters (Percentage of samples with PSI value ≥75%, standard deviation of PSI value >0.1 and average PSI value > 0.05).

Every AS event was assigned a unique identifier by combining the gene symbol, the ID number in the SpliceSeq database and splicing type. For instance, in the identifier term “CSAD-21952-ES”, the gene symbol is CSAD, ID number is 21,952 and splicing type is ES.

Identification of overall survival related AS events, functional enrichment analysis and gene network construction

We identified overall survival (OS) related AS events using univariate Cox regression analysis with a P value <0.05. Upset plot which is similar to a Venn diagram, was introduced to depict the intersections between the seven types of AS events in HCC (Khan & Mathelier, 2017). The corresponding genes of OS-related AS events were sent for functional enrichment analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways with P value <0.05 and false discovery rate (FDR) <0.05 were considered significantly. The top significant pathways in KEGG and GO were showed with bar plots. Gene network analysis was performed by inputting the corresponding genes of 800 most significant OS-related AS events to Cytoscape and Reactome FI (version 3.7.1), further, hub genes were selected at the same time.

Construction of the prediction model based on AS events

Least absolute shrinkage and selection operator (LASSO) regression is suitable for the reduction of high-dimensional data to avoid model overfitting (Sauerbrei, Royston & Binder, 2007). When we performed LASSO regression, some variables were eliminated through penalty rules and potential predictors with non-zero coefficients were finally leaved (Gao, Kwan & Shi, 2010). We determined the penalty parameter lambda by the cross-validation using the glmnet package. The optimal lambda value corresponding to the minimum value of the cross-validation error mean was identified to determine the potential OS-related AS events (Tibshirani, 1997). We selected the optimal OS-related AS events in each type respectively with nonzero coefficients in the LASSO regression, and constructed the prediction models by using multivariate Cox regression analysis. The final prediction model (all types) was constructed by combing seven types of AS events also selected via LASSO regression.

The prediction model was established according to risk score, and the risk score was calculated by the PSI value of each AS events and the corresponding regression coefficient (lnHR) generated from multivariate Cox regression analysis. The formula is as follows: $Riskscore={\sum }_i^{n}PSIi\ast \beta i$, where β is the regression coefficient.

Prognostic value evaluation of the prediction model

The HCC patients were divided into low and high risk groups according to the median value of the risk score and Kaplan–Meier survival analysis was performed to compare the OS rate between the two groups within five years. The p-values were computed using log-rank test. Furthermore, ROC curves of 1, 3, 5 years were generated to compare the predictive accuracy of each prediction model.

Construction of the regulation network between SFs and AS events

SF data was obtained from the SpliceAid2 database (http://www.introni.it/splicing.html). The association of SF expression with OS was analyzed using univariate Cox regression. The relationship between the PSI value of OS-related AS events and SF expression were analyzed using Spearman correlation analysis. Regulation network was generated by Cytoscape and Reactome FI (version 3.7.1).

Statistical analysis

The R software 3.5.0 was utilized for all statistical analysis. A P value <0.05 was deemed to be of statistical significance.

Overview of AS events in HCC cohort

A flowchart of our study design is showed in Fig. 1. We detected a total of 34,163 AS events in 8985 genes in 335 HCC patients, comprised of 12,327 ES events in 5,343 genes, 8087 AT events in 3,532 genes, 6,352 AP events in 2,566 genes, 2,666 AA events in 1,937 genes, 2,331 AD events in 1,663 genes, 2,263 RI events in 1,561 genes and 137 ME events in 135 genes (Table 1). We can find that one single gene undergoes up to six types of AS events from the UpSet plot (Fig. 2A). In addition, ES was the most common among seven types of AS events.

The OS-related AS events of HCC

A total of 3,482 OS-related AS events in 2,203 genes were detected using univariate Cox analysis, including 1,264 ES events in 988 genes, 887 AT events in 534 genes, 643 AP events in 411 genes, 217 AA events in 203 genes, 231 AD events in 205 genes, 224 RI events in 196 genes and 16 ME events in 16 genes (Table 1). A volcano plot of these AS events was provided in Fig. 2B. With the display of UpSet plot, one single gene could have up to four types OS-related AS events (Fig. 2C). The top 20 significant OS-related AS events (if available) for each AS type were showed by bubble plots in Figs. 3A–3G. Obviously, most of AS events in ES were favorable prognostic elements. However, most of these AS events in RI, AA and AD were adverse prognostic elements.

We established a gene interaction network by sending the corresponding genes of the top 800 most significant OS-related AS events to Cytoscape. The results indicated that UBB, UBE2D3, SF3A1 were the hub genes of this gene network (Fig. 4A). Furthermore, 2,203 genes from 3,482 OS-related AS events were used for KEGG and GO enrichment analysis to explore the pathways and biological functions of the OS-related AS genes. The top 10 significant terms of GO enrichment analysis were presented in Fig. 4B. Such as spliceosomal complex, adherens junction in cellular component (CC); damaged DNA binding, cell adhesion molecule binding in molecular function (MF); protein targeting, actin cytoskeleton reorganization in biological process (BP). A total of nine KEGG pathways were identified, such as base excision repair, pyruvate metabolism and PPAR signaling pathway (Fig. 4C).

Construction of the prediction model for HCC patients

We used LASSO regression to select the top 10 optimal OS-related AS events (if available) and then construct the prediction model (Fig. 5). Three OS-related AS events were selected for ES; four OS-related AS events for ME; five OS-related AS events for AP and AT; six OS-related AS events for AA, AD and RI; nine OS-related AS events for the final prediction model (Table 2). Risk scores were computed according to the selected AS events, and HCC patients were divided into low and high risk groups on the basis of the median value of risk scores. The distribution of survival status in low and high risk groups, risk score curves and the PSI value heat map of the AS events for eight prediction models were visualized in Fig. 6. The results of Kaplan–Meier survival analysis showed that all of the eight prediction models possess significant ability to predict the prognosis of HCC patients between low and high risk group (Fig. 7). However, according the results of ROC curves, the final prediction model exhibited the most powerful predictive efficiency than other models established by one single AS type with the maximum AUCs of 0.878, 0.843, 0.821 in 1, 3, 5 years ROC curves (Fig. 8).

Moreover, we assessed the prognostic value of risk score and other clinicopathological characteristics using univariate and multivariate Cox regression analysis. Univariate analysis suggested that advanced clinical stage (HR = 2.053, 95% CI [1.553–2.714], p < 0.001), high T classification (HR = 1.962, 95% CI [1.515–2.542], p < 0.001) and high risk score (All) (HR = 1.105, 95% CI [1.076–1.135], p < 0.001) were associated with poor prognosis (Fig. 9A). Upon multivariate analysis, risk score (All) (HR = 1.101, 95% CI [1.069–1.135], p < 0.001) served as the risk factor to independently predict OS for HCC patients (Fig. 9B).

The regulation network of SFs and OS-related AS events

SF is the executor of AS events and aberrant expression of SF was related to oncogenesis process. A total of 71 SFs data were collected from the SpliceAid2 database (http://www.introni.it/splicing.html). We identified 21 SFs associated with OS of HCC patients using univariate analysis. Then Spearman correlation analysis was performed to determine the correlation between the PSI value of OS-related AS events and SF expression. The significant correlations (—R—≥0.6, p < 0.001) were selected to construct the regulation network (Fig. 10A). The regulation network consists of 29 OS-related AS events, of which 16 were adverse AS events (red dots) and 13 were favorable AS events (green dots), were significantly correlated with the 13 SFs (blue dots). We can find that all of the SFs were correlated with multiple AS events and played opposite roles in regulating different AS events. Similarly, a part of the AS events could be regulated by different SFs. This phenomenon partly explains that the same transcript can produce multiple different splicing events. In addition, we can detect that the adverse AS events were positively correlated with SF expression (red lines), whereas the favorable AS events were negatively correlated with SF expression (green lines).

The top four most significant correlations between SFs and OS-related AS events (—R— ≥0.7) are shown in Figs. 10D, 10E, 10H and 10I. The top two counterpart SFs were HNRNPL and HNRNPH3, and the expression of HNRNPL and HNRNPH3 in tumor sample were significantly higher than adjacent normal sample (Figs. 10B and 10F). HCC patients were classified into low and high risk groups according to the median value of the two SFs expression. Also, Kaplan–Meier survival analysis showed statistical difference in survival time between the two groups (Figs. 10C and 10G).