Prognostic index of aberrant mRNA splicing profiling acts as a predictive indicator for hepatocellular carcinoma based on TCGA SpliceSeq data.

Alternative splicing in tumor cells may be used as a molecular marker for the differential diagnosis of certain tumor types and assessment of prognosis. The aim of the present study was to investigate the associations among alternative splicing events, splicing factors, and the survival of patients with hepatocellular carcinoma (HCC). The alternative splicing event profiles of 371 patients with HCC were downloaded from The Cancer Genome Atlas (TCGA) SpliceSeq data, and the percent‑splice‑in value for each splicing event was calculated. The association between alternative splicing events and overall survival was evaluated. The most significant prognosis‑related splicing events were used to build up a prognostic index (PI). A total of 3,082 survival‑associated alternative splicing events were detected in HCC. The final PI based on all of the most significant candidate alternative splicing events exhibited better performance in distinguishing good or poor survival in patients compared to the PI based on a single type of splicing event. Receiver operating characteristic curves confirmed the high efficiency of the PI in predicting the survival of HCC patients, with an area under the curve of 0.914. The overexpression of 32 prognosis‑related splicing factor genes could also predict poor prognosis in patients with HCC. In conclusion, the constructed computational prognostic model based on HCC‑specific alternative splicing events may be used as a molecular marker for the prognosis of HCC.

Introduction

Alternative splicing, also referred to as differential splicing, refers to the process of generating various isoforms of an mRNA precursor by several approaches, i.e., by selecting diverse splice site combinations. The resulting protein products may exhibit distinguishable or antagonistic functional and structural characteristics, or they may yield various phenotypes in the same cell type due to differences in their expression levels (1,2).

The splicing isoforms of certain genes may serve as drivers for the development of tumors (3-7). Additionally, the genomes of tumor cells have more extensive alternative splicing mechanisms compared with normal cells. In most types of tumors, the expression products of genes may vary due to changes in alternative splicing, which is associated with tumor development, metastasis (8-12) and resistance to treatment (7,13-15). Additionally, alternative splicing events in tumors may act as diagnostic and prognostic molecular biomarkers (16-19), and they may also be used as targets for cancer treatment (20-24).

Alternative splicing events are regulated by splicing factors that act as oncogenes or tumor suppressor genes (25). The development of drugs targeting splicing factors is a new research focus in cancer therapy (26-32), and several studies have attempted to distinguish cancers based on the presence of abnormal RNA splicing (10,33). Furthermore, disturbances of alterative splicing are reportedly involved in the development of cancers. In addition, several survival-associated splicing events that are unique to cancer have been identified (34). RNA-binding proteins, which function as regulatory elements of splicing events, have also been associated with cancer development (35-37).

Hepatocellular carcinoma (HCC) is the fifth most common type of cancer with the third highest mortality rate worldwide. Its occurrence has been associated with viral infection, aflatoxin exposure, and other factors; in addition, a number of studies have demonstrated that alternative splicing also plays an important role in HCC (38-43). However, these studies focused on the associations between the alternative splicing events of only a few genes and HCC, whereas only few studies have analyzed alternative splicing events in HCC using large datasets. Tremblay et al (44) analyzed RNA sequencing (RNA-Seq) data from The Cancer Genome Atlas (TCGA) and found a large number of differential alternative splicing events between HCC and corresponding paracancerous tissues. Hepatitis B and C virus infections also affect alternative splicing events in HCC (44). Accordingly, alternative splicing is considered to play an important role in the occurrence of HCC. However, to the best of our knowledge, no study has yet examined the association between mRNA splicing and the prognosis of HCC based on TCGA SpliceSeq data. Therefore, by using TCGA SpliceSeq data, the present study analyzed the associations among alternative splicing events, splicing factors, and the survival of patients with HCC, with the aim of identifying splicing events that may serve as new molecular targets for the prognosis of HCC.

Materials and methods

Assortment of alternative splicing event data

The alternative splicing event profiles of HCC patients were downloaded from TCGA SpliceSeq (45), which is a resource for the investigation of transcript splicing patterns and splicing event details based on TCGA covering quantified introns or exons. Information on percent-splice-in (PSI), the ratio of normalized read counts indicating the inclusion of a transcript element over the total normalized reads for that event, was collected from the database. PSI was calculated as the ratio of reading densities of inclusions to the sum of the reading densities of inclusion and exclusion. PSI values range from 0 to 100%. Only samples with PSI values >90% were downloaded. Simultaneously, clinical data were also obtained from TCGA. A total of 7 different alternative splicing events were obtained: Exon skips (ESs), retained introns (RIs), mutually exclusive exons (MEs), alternate donor sites (ADs), alternate acceptor sites (AAs), alternate promoters (APs), and alternate terminators(ATs).

Survival analysis and production of a prognostic signature

A total of 371 HCC patients in TCGA were used to select survival-related alternative splicing events. The association between alternative splicing events and overall survival (OS) was evaluated using univariate Cox regression analysis. The most significant prognostic alternative splicing events (P<0.0001) were subjected to multivariate Cox regression analysis. The area under the curve (AUC) of a time-dependent receiver operating characteristic (ROC) curve was also calculated, which has been used widely to compare the ability of prognostic predictors. All survival analyses were conducted using the 'survival' and 'survival ROC' packages in R software. The OS rates of patients with HCC grouped by high- and low-risk alternative splicing events were plotted using Kaplan-Meier plots and the differences between groups were analyzed using log-rank tests.

UpSet plot and gene-annotation enrichment analysis

An UpSet plot, a novel visualization technique for the quantitative analysis of interactive sets, was used to analyze the intersections between the 7 types of alternative splicing. Gene-annotation enrichment analyses were used in the 'clusterProfiler' package in R to annotate and visualize biological process terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of survival-associated alternative splicing genes based on the criterion of an adjusted P-value of <0.05.

Splicing factor genes

Splicing factor genes were collected from the SpliceAid 2 database (http://193.206.120.249/splicing_tissue.html). Subsequently, the expression profiles of splicing factor genes were extracted from TCGA. Cox univariate regression analysis was performed, and genes significantly associated with OS (P<0.001) were retained. An analysis of the 32-gene splicing factors was also incorporated by using the ProgGene (http://genomics.jefferson.edu/proggene/) database (46). The database helps users perform a combined analysis for a list of genes and generate a prognostic signature based on the imputed genes using Cox proportional hazard analysis. The 32 splicing factors were input and the TCGA database was selected for survival analysis. Student t-test analysis was conducted to estimate the differences of splicing factors between tumor and non-tumor samples.

Splicing correlation network construction and hub splicing factor genes

The correlation between the expression of prognosis-related splicing factor genes and the PSI values of the most significant survival-associated alternative splicing events were analyzed by Pearson's tests. The splicing events-factors axes were submitted to correlation network construction for P-values <0.05 and a Pearson's correlation coefficient >0.3. The correlation plots were generated by Cytoscape (version 3.5.1). Hub splicing factor genes were identified according to the numerical degree of each node and edge. The expression level of hub splicing factors between tumor and non-tumor tissues were calculated based on data from TCGA. The data were transformed to the form of log2(count+1) and visually displayed by OriginPro 2017 (OriginLab Corp.). Furthermore, calibration plots were generated to evaluate model calibration by using 'rms' package in R software (https://CRAN.R-project.org/package=rms).

Results

Identification of mRNA splicing event profiles in HCC

The alternative splicing event profiles of 371 HCC patients were obtained from TCGA SpliceSeq. In total, 2,6210 mRNA splicing events in 7,727 genes were collected, which comprised 8,524 ESs in 4,120 genes, 1,306 RIs in 1,879 genes, 4,359 APs in 1,775 genes, 7,736 ATs in 3,379 genes, 1,675 ADs in 1,205genes, 1,936 AAs in 1,451 genes, and 102 MEs in 101 genes. These results also indicated that one gene may have several types of mRNA splicing events, and ES was the predominant type.

Prognosis-related alternative splicing events in HCC

The association of these alternative splicing events with OS was studied by univariate survival analysis, and a total of 3,082 survival-associated alternative splicing events were detected in HCC (P<0.05). We observed that one gene may have two or more events that are significantly associated with survival. Thus, an intersection visualization plot was generated via UpSet (Fig. 1). The top significant survival-associated genes in HCC (P<0.0001) were piped to annotate and visualize biological process terms and KEGG pathways. The biological process terms of these genes were most enriched in small molecule catabolic process, organic acid catabolic process, and viral gene expression (Fig. 2). KEGG pathway analysis indicated that these genes were mainly enriched in ribosome, purine metabolism, protein processing in endoplasmic reticulum, and autophagy.

Figure 1

UpSet plot of survival-related alternative splicing events. Venn diagram of the overlap between the 7 types of survival-associated alternative splicing events in HCC. HCC, hepatocellular carcinoma; AA, alternate acceptor sites; AD, alternate donor sites; AP, alternate promoters; AT, alternate terminators; ES, exon skipping; ME, mutually exclusive exons; RI, retained introns.

Figure 2

Pathway analysis of survival-related alternative splicing genes. (A) Biological process analysis of survival-related alternative splicing genes. (B) Significantly enriched annotation from KEGG pathway analysis of survival-related alternative splicing genes. KEGG, Kyoto Encyclopedia of Genes and Genomes; ER, endoplasmic reticulum.

Prognostic indicators for HCC patients

The top significant survival-associated alternative splicing events in the 7 types of alternative splicing were selected as candidate prognostic factors. Multivariate Cox regression with prognostic model construction was then applied to develop a computational model for HCC prognosis based on each individual splicing type. HCC patients were divided into subgroups with significantly different prognoses (Fig. 3) based on the 7 types of alternative splicing (Table I). These independent prognostic alternative splicing events were analyzed further to build the final prognostic computational model using a PI for HCC. Time-dependent ROC curve analyses were also applied to compare the efficiency of these prognostic models (Figs. 4 and 5). The final PI containing the most significant candidate alternative splicing events exhibited the best performance at distinguishing favorable or poor survival in patients (Fig. 6). Patients in the high-risk group had significantly worse survival compared with those in the low-risk group [hazard ratio (HR)=12.904, 95% confidence interval (CI): 7.513-22.161, P<0.001; Fig. 6A and B]. Furthermore, the PI remained an independent prognostic indicator for HCC patients in multivariate analyses after other clinicopathological characteristics were adjusted (HR=16.541, 95% CI: 7.660-35.717, P<0.001; Table II). ROC curves confirmed that the final prognostic predictor with all types of alternative splicing had the highest efficiency (Fig. 6C) with an AUC of 0.914. The specific genes involved in the final models are listed in Table III. The calibration curves exhibited good consistency between the predicted and actual survival probability (Fig. 6D).

Figure 3

Kaplan-Meier analyses indicate that HCC patients can be stratified into low- risk and high-risk groups based on the 7 types of alternative splicing. (A) AA; (B) AD; (C) AP; (D) AT; (E) ES; (F) ME; and (G) RI. HCC, hepatocellular carcinoma; AA, alternate acceptor sites; AD, alternate donor sites; AP, alternate promoters; AT, alternate terminators; ES, exon skipping; ME, mutually exclusive exons; RI, retained introns.

Table I

Prognostic signatures based on each type of splicing event.

Type of splicing event	Algorithm	HR (95% CI)	ROC curve
Alternate acceptor sites (AA)	'SWI5-87732-AA' * 0.03740 + 'FDPS-8074-AA' * 0.02940 + 'COMT-61101-AA' * 0.31109	2.334 (1.644-3.314)	0.743
Alternate donor sites (AD)	'VPS28-85610-AD' * 1.70e-01 + 'UBE2D3-70148-AD' * (−1.57e-02) + 'TXNDC17-38768-AD' * 9.27e-02 + 'TPT1-25803-AD' * 1.62e+01 + 'SUMF2-79789-AD' * 3.92e-02 + 'SMIM19-83740-AD' * 2.69e-02 + 'RPL26-39180-AD' * 6.06e-01 + 'NDUFB11-88904-AD' * 4.87e-02 + 'IRF3-50997-AD' * (−8.62e-03) + 'F10-26354-AD' * 2.21e-02 + 'CCDC90B-18085-AD' * 2.02e-02 + 'C1orf159-13-AD' * (−1.12e-02) + 'ADCK4-49949-AD' (-4.56e-02)	2.423 (1.700-3.455)	0.718
Alternate promoters (AP)	'RPL27-41172-AP' * 0.21696 + 'RCL1-85779-AP' * −0.00550 + 'NUDT6-70523-AP' * 0.01129 + 'MXI1-13080-AP' * 0.00857 + 'DUSP22-75126-AP' * 0.07876 + 'CYB5A-45803-AP' * 0.49517 + 'CHPF-57719-AP' * 0.07805	2.687 (1.891-3.818)	0.728
Alternate terminators (AT)	'ZSCAN20-1679-AT' * 0.01077 + 'ZNF706-84737-AT' * 0.06539 +'WWOX-37672-AT' * −0.03181 + 'NUCB2-14523-AT' * −0.07215 +'NEIL3-71227-AT' * 0.01306 + 'MORN1-252-AT' * 0.00837 + 'FRMD4A-10810-AT' * 4.59555 + 'CRISPLD2-37867-AT' * 0.24507 + 'ACAT1-18599-AT' * 0.69847	2.850 (2.001-4.057)	0.709
Exon skips (ES)	'USF2-49098-ES' * 0.03030 + 'TRAPPC6A-50409-ES' −0.03711 + 'TMEM150A-54305-ES' −0.01822 + 'TMEM120B-24894-ES' * 0.00897 + 'TCF3-46542-ES' * −0.02326 + 'SSR4-90499-ES' * 0.95791 + 'SMIM7-48191-ES' * 0.26393 + 'SHBG-39015-ES' * −0.02339 'PTRH1-87649-ES' * −0.64609 + 'PPP6R2-62825-ES' * −0.03581 + 'PPP2CB-83304-ES' * −3.17592 + 'PGAP3-40670-ES' * 0.06043 + 'PCYT2-44234-ES' * 0.12079 + 'OCEL1-48244-ES' * −0.04796 + 'SOC2-52105-ES' * 0.01667 + 'IRF3-50995-ES' * −0.03249 + 'IRF3-50990-ES' * 0.04115 + 'FKBP8-48448-ES' * −0.71033 + 'ECHDC2-3037-ES' * 0.04796 + 'CYTH1-43890-ES' * 0.01764 + 'CHRD-67959-ES' * −0.01020 + 'C16orf13-32924-ES' * 0.03301 + 'C14orf2-29536-ES' * 0.38460 + 'BAX-50838-ES' * −0.35118 + 'ARHGEF10L-862-ES' * 0.01344 + 'AFMID-94694-ES' * −0.01010 + 'ACAA1-64027-ES' * −0.71363	4.938 (3.439-7.090)	0.898
Mutually exclusive exons (ME)	'SLC39A14-140283-ME' * (−0.00563) + 'H2AFY-96931-ME' * 0.02881	1.878 (1.320-2.672)	0.641
Retained introns (RI)	'RPL13-38091-RI' * 0.62076 + 'ROMO1-59223-RI` * 0.07039 + 'NUDT22-16590-RI' * 0.01621 + 'NAA60-33527-RI' * (−0.00942) + 'MRPL52-26642-RI' * 0.01287 + 'MDK-15571-RI' * 0.09757 + 'MARVELD2-72372-RI' * (−2.62393)	2.895 (2.043-4.102)	0.745

HR, hazard ratio; CI, confidence interval; ROC, receiver operating characteristic.

Figure 4

ROC curves of the PI built using the 7 types of alternative splicing in HCC. (A) AA; (B) AD; (C) AP; (D) AT; (E) ES; (F) ME; and (G) RI. ROC, receiver operating characteristic; AUC, area under the ROC curve; PI, prognostic index; HCC, hepatocellular carcinoma; AA, alternate acceptor sites; AD, alternate donor sites; AP, alternate promoters; AT, alternate terminators; ES, exon skipping; ME, mutually exclusive exons; RI, retained introns.

Figure 5

The 7 prognostic models with single prognosis-relevant splicing events. (A) Alternate acceptor sites; (B) alternate donor sites; (C) alternate promoters; (D) alternate terminators; (E) exon skipping; (F) mutually exclusive exons; and (G) retained introns.

Figure 6

Construction of the PI of HCC patients based on all alternative splicing events. (A) Kaplan-Meier plot indicating that patients in the high-risk group had significantly shorter OS compared with those in the low-risk group. (B) OS status and survival duration of HCC patients. (C) Predictive value of the PI for clinical outcome by ROC curve. (D) Calibration curves for predictions of OS. ROC, receiver operating characteristic; AUC, area under the ROC curve; PI, prognostic index; HCC, hepatocellular carcinoma; OS, overall survival.

Table II

Univariate and multivariate analyses of OS in HCC patients from TCGA by Cox regression analysis.

	Univariate analysis		Multivariate analysis
Variables	HR (95% CI)	P-value	HR (95% CI)	P-value
Age (≥60/<60 years)	1.212 (0.854-1.720)	0.281
Sex (male/female)	0.817 (0.573-1.164)	0.262
Pathological stage (III-IV/I-II)	2.446 (1.687-3.545)	<0.001	1.018 (0.138-7.527)	0.986
T stage (T3-T4/T1-T2)	2.537 (1.783-3.609)	<0.001	1.757 (0.238-12.981)	0.581
Lymph node metastasis (yes/no)	1.999 (0.490-8.161)	0.334
Distant metastasis (yes/no)	4.033 (1.267-12.834)	0.018	0.755 (0.179-3.177)	0.701
Histological grade (G3-G4/G1-G2)	1.119 (0.780-1.604)	0.542
Tumor status during follow-up (tumor/tumor-free)	2.366 (1.623-3.447)	<0.001	1.567 (0.979-2.508)	0.061
Vascular invasion (micro+macro/none)	1.351 (0.892-2.047)	0.155
Relative family cancer history (yes/no)	1.181 (0.818-1.704)	0.375
Prognostic index (high/low)	12.904 (7.513-22.161)	<0.001	16.541 (7.660-35.717)	<0.001

OS, overall survival; HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; HR, hazard ratio; CI, confidence interval.

Table III

Prognostic index for HCC patients based on alternative splicing events.

Symbol	AS ID	Splice type	Coefficient	z score	P-value
ZFAND6	32171	AA	−8.51E-01	−2.83	0.0047
SWI5	87732	AA	8.52E-01	2.46	0.01392
BIRC2	18447	AA	−1.28E+01	−4.66	3.20E-06
RPS14	74096	AD	−2.55E+00	−4.38	1.20E-05
PGRMC2	70577	AD	−3.42E+01	−4.58	4.70E-06
CCDC90B	18085	AD	8.30E-01	2.57	0.0102
APOC1	50361	AD	4.11E+00	3.32	0.00089
RTN4	53584	AP	−4.26E-01	−3.29	0.001
RPL23A	39967	AP	−9.34E+01	−2.57	0.01003
RPL23A	39965	AP	−3.06E+00	−3.44	0.00059
RCL1	85779	AP	4.21E-01	2.92	0.00347
NUDT6	70523	AP	2.19E-01	2.16	0.03054
NUDT6	70521	AP	6.16E-01	3.07	0.00216
MXI1	13080	AP	1.53E-01	1.49	0.13594
MXI1	13079	AP	4.76E-01	2.67	0.00752
IAH1	52628	AP	1.41E+01	4.98	6.30E-07
DUSP22	75126	AP	−2.38E-01	−1.76	0.07903
DUSP22	75125	AP	9.50E-01	3.53	0.00041
CYB5A	45803	AP	1.29E+00	3.29	0.00099
CHPF	57718	AP	−1.64E+00	−2.21	0.02698
CARKD	26251	AP	−2.49E-01	−1.71	0.08786
IAH1	52629	AP	4.82E-01	3	0.00273
ZSCAN20	1679	AT	−7.88E-01	−3.12	0.00182
ZNF706	84738	AT	−1.42E+01	−5.76	8.40E-09
WWOX	37672	AT	−5.47E+00	−7.13	9.70E-13
TROAP	21551	AT	3.63E-01	1.56	0.11927
RAPH1	57075	AT	1.48E+00	1.83	0.0673
RAPH1	57074	AT	5.54E-01	2.8	0.00504
NUCB2	14523	AT	−9.05E+00	−2.22	0.0262
FRMD4A	10810	AT	5.09E+00	5.4	6.50E-08
CRISPLD2	37867	AT	8.61E-01	3.72	0.0002
CLSPN	1731	AT	−8.03E-01	−1.74	0.08206
BCAM	50347	AT	−3.64E-01	−1.6	0.11049
ACAT1	18599	AT	2.88E+00	4.6	4.30E-06
XPC	63523	ES	2.96E-01	1.79	0.0731
WIPI2	78656	ES	−3.47E-01	−1.98	0.04758
TRAPPC6A	50409	ES	−3.27E+00	−4.69	2.70E-06
TPRA1	66613	ES	−4.71E+00	−1.84	0.06575
TMEM150A	54305	ES	1.66E+00	2.85	0.00443
TMEM141	88206	ES	−3.62E-01	−3.7	0.00022
TCF3	46542	ES	−1.57E+00	−4.66	3.10E-06
STRA13	44266	ES	2.97E+00	4.89	1.00E-06
SSR4	90499	ES	2.39E+02	6.3	3.00E-10
SLC7A9	48907	ES	3.61E-01	3.17	0.0015
SLC25A45	16828	ES	3.71E-01	2.05	0.04004
SLAIN2	69214	ES	−6.78E-01	−2.89	0.00385
SHBG	39015	ES	−8.49E-01	−7.07	1.50E-12
PTRH1	87649	ES	−6.58E+01	−5.47	4.40E-08
PPP2CB	83304	ES	−5.52E+02	−5.72	1.10E-08
PGAP3	40670	ES	2.10E+01	3.32	0.00089
OCEL1	48244	ES	−5.74E+00	−2.97	0.00298
MTA1	29647	ES	−3.93E-01	−2.25	0.02464
MPPE1	44648	ES	−3.09E-01	−2.34	0.01921
MPND	46796	ES	−6.87E+00	−2.12	0.034
LTBP3	16865	ES	−4.93E-01	−1.87	0.06167
JOSD2	51205	ES	−5.21E+00	−2.51	0.0119
ISOC2	52106	ES	−5.98E+00	−3.29	0.00099
IRF3	50991	ES	−7.04E+01	−2.94	0.00327
IRF3	50990	ES	5.08E-01	2.41	0.01599
GUK1	10185	ES	−2.01E+01	−3.77	0.00016
ECHDC2	3037	ES	6.27E-01	3.96	7.50E-05
DGKZ	15546	ES	−8.54E-01	−2.5	0.01235
CSAD	21952	ES	8.72E-01	3.91	9.30E-05
CLASRP	50393	ES	−1.15E+00	−2.99	0.00281
CHRD	67959	ES	−5.44E-01	−3.86	0.00011
C17orf49	38825	ES	2.01E+01	3.52	0.00043
C16orf13	32921	ES	2.25E+00	2.39	0.01699
C16orf13	32919	ES	−3.45E+00	−3.98	6.80E-05
BAX	50838	ES	−4.99E+01	−3.86	0.00011
ARHGEF10L	862	ES	4.31E-01	1.83	0.06684
APOC1	99361	ES	−2.72E-01	−3.56	0.00037
AFMID	94694	ES	2.79E-01	1.98	0.04811
AFMID	94690	ES	4.05E-01	2.62	0.00888
AFMID	43800	ES	−1.21E+00	−4.21	2.50E-05
AFMID	43798	ES	7.48E-01	3.04	0.00236
AFMID	43795	ES	−5.04E-01	−4.01	6.00E-05
ACY1	65150	ES	3.56E+01	1.76	0.07779
SLC39A14	140283	ME	3.82E-01	2.9	0.00372
H2AFY	96931	ME	2.26E+00	3.09	0.00198
CYP4F3	48101	ME	−3.19E-01	−2.94	0.00327
VPS28	85600	RI	−7.11E-01	−2.25	0.02443
ROMO1	59223	RI	7.85E-01	2.64	0.00822
NUDT22	16590	RI	2.71E+00	5.83	5.70E-09
MRPL52	26642	RI	1.11E+00	4.03	5.60E-05

AS, alternative splicing; HCC, hepatocellular carcinoma; AA, alternate acceptor sites; AD, alternate donor sites; AP, alternate promoters; AT, alternate terminators; ES, exon skipping; ME, mutually exclusive exons; RI, retained introns.

Survival-associated splicing factor genes

A total of 71 splicing factor genes were identified, and the associations between the expression profiles of these genes and OS were assessed based on the data obtained from TCGA, resulting in the identification of 32 prognosis-related splicing factor genes. Interestingly, overexpression of all these 32 genes may represent a risk factor for the survival of patients with HCC, as all HRs were >1 (P<0.05) (Fig. 7A). We found that the integration of the 32 SFs could also provide an effective survival risk stratification for HCC patients (Fig. 7B). HCC patients with high-risk scores for the splicing factors-based classifier had worse OS compared with those who had low-risk scores. Next, we characterized the splicing-regulatory network of the 152 (P<0.0001) most significant survival-associated alternative splicing events and 32 prognosis-related splicing factor genes in HCC (Fig. 8). Additionally, 14 core interacting genes were identified in the splicing-regulated network: SRSF7, SRSF10, SRSF1, RBMX, HNRNPC, YBX1, PTBP1, HNRNPA3, HNRNPU, HNRNPL, SRSF9, SRSF2, DAZAP1 and SFPQ (Fig. 9). Furthermore, these genes were found to not only possess prognostic significance, but to also have indispensable functions in the initiation of HCC (Figs. 10 and 11).

Figure 7

Prognostic value of survival-associated splicing factor genes. (A) Forest plots of the HRs of survival-associated splicing factor genes in HCC. (B) The integration of the 32 splicing factor genes may achieve effective survival prediction. HR, hazard ratio; CI, confidence interval; HCC, hepatocellular carcinoma.

Figure 8

Survival-associated splicing factors and splicing correlation network in hepatocellular carcinoma. The red circles indicate risk alternative splicing events and the blue circles indicate protective alternative splicing events. The blue circles represent splicing factor genes. The red and blue lines represent positive and negative correlations, respectively.

Figure 9

Bar chart of the top 30 splicing factor genes in hepatocellular carcinoma.

Figure 10

Heatmap of hub splicing factor gene expression levels in HCC and adjacent non-tumor tissues. HCC, hepatocellular carcinoma.

Figure 11

Expression patterns of 14 hub splicing factor genes in HCC and paired non-tumor samples. Each red dot represents a distinct tumor sample and each blue dot represents a non-tumor sample. HCC, hepatocellular carcinoma. *P<0.05.

Discussion

In the present study, the PSI value and clinical outcome of patients with HCC were integrated using a novel methodology. Furthermore, prognostic signatures were developed based on splicing events to promote the development of precision treatment. To determine the underlying mechanism, we investigated the potential regulatory associations between prognosis-related splicing factors and splicing events, which provided deeper insights into the phenotype of splicing events in HCC.

With the advantage of high-throughput RNA-Seq, TCGA provides multiple sources for the investigation of whole-genome or transcriptome analyses, including the exploration of genome splicing events. Previously, Xue et al (47) and Zhu et al (48) conducted SpliceSeq analyses using TCGA (49) to generate alterative splicing profiles for non-small-cell lung cancer and ovarian cancer patients, in order to develop prognostic models with several splicing events. SpliceSeq is a Java program providing a clear view of the inclusion level of each exon and splice junction. Recently, Ryan et al (45) extended the methodology of SpliceSeq and calculated the PSI for each potential splicing event across 33 types of cancer to establish TCGA SpliceSeq database; the PSI data for HCC generated in the present study were obtained from this database. Hence, to the best of our knowledge, this is the first group to integrate PSI values from SpliceSeq with TCGA clinical prognostic parameters for the comprehensive investigation of the prognostic value of alternative splicing events in HCC. This algorithm will facilitate the discovery of novel splicing events that may be used for the prognosis of HCC and provide insights into the role of splicing events at the genome-wide level.

On the basis of the analysis of the tumor tissue data of 371 HCC cases in TCGA SpliceSeq, a large number of alternative splicing events were identified. These events do not only affect protein diversity by altering the amino acid sequences of the translated products (e.g., ES can delete a portion of the coding sequence), but they can also affect protein function by altering the yield of the translated products (e.g., AP can change the efficiency of translation). For example, the functions of the common clinical tumor suppressor genes TP53, ARID1A, PTEN and PI3K, and the proto-oncogenes MET and NOTCH1, can be changed by alternative splicing events (21). The common apoptosis-related genes Bcl-X and MLC1 can also express different products due to the occurrence of alternative splicing, showing two distinct functions of promoting apoptosis and inhibiting apoptosis (3,4,50,51). In the currently popular molecular targeting therapy for tumors, the target sequence can be lost from cancer cells due to alternative splicing, leading to drug resistance (7,13,21). In the present study, ES was the most common of the 7 alternative splicing events in HCC. Previous studies have demonstrated that the frequency of alternative splicing events is not the same in different cancers, but the frequency of ES is generally the highest (52). Although the mechanism has not yet been fully elucidated, evidence has shown that the occurrence of splicing events is not random. Chen et al (53) analyzed RNA-Seq data from 30 types of tumors and found that the frequency of alternative splicing events is markedly higher in tumors compared with that in corresponding normal tissues. These alternative splicing events are often accompanied by premature termination codons, and the probability of a premature termination codon due to alternative splicing is markedly higher in a tumor suppressor gene rather than in a proto-oncogene in cancer tissue. In the gene functional enrichment analysis, it was observed that genes with survival-associated splicing events were significant correlated with 'small molecule catabolic processes', 'organic acid processes' and 'viral gene expression'. The gene functional enrichment analysis results revealed that these genes may be involved in the metabolism of HCC. The liver is the metabolic center of the body. Hence, the onset and progression of HCC is also inseparable from metabolic abnormalities. Previously, Tremblay et al reported unique alternative splicing patterns in HBV-associated HCC, HCV-associated HCC, HBV&HCV-associated HCC, and virus-free HCC based on TCGA dataset. They also found the signatures of genes for which AS is dysregulated in different types of HCC (44). These findings also validated that the process of splicing is closely associated with viral infection status.

Accordingly, alternative splicing events may be driving the occurrence and progression of cancer. Therefore, the Cox regression method was used to analyze the association between alternative splicing events and the survival of HCC patients. The predictive power of alternative splicing on the prognosis of HCC was analyzed by considering the ROC curve. The results demonstrated that several alternative splicing events were associated with the survival of HCC patients. Thus, we attempted for the first time to use multiple alternative splicing events for the prognosis of HCC patients and obtained satisfactory results. To date, a number of studies have demonstrated that an alternative splicing event alone may be used as a biomarker for tumor diagnosis and prognosis. However, due to the heterogeneity of tumors, using a single molecule as a marker is often insufficient. Therefore, the combination of multiple molecular events for the prognosis of HCC may improve the sensitivity of diagnosis and prognosis. This follows a trend in current clinical research. Similarly, previous studies have reported several prognostic signatures based on different types of molecular events. For example, Qiao et al proposed an eight-gene signature for HCC patients' survival prediction (54). However, the AUC of their model only reached 0.77. Liao et al also developed an effective prognostic index (PI) for predicting the survival of HCC patients based on microRNAs (55). However, the AUC also only reached 0.687. Hence, our PI exhibits better discriminatory power. In the present study, the discriminatory power and model calibration were estimated to comprehensively evaluate the model accuracy. More importantly, the prognostic signature based on splicing events may provide novel insight into the monitoring of HCC patients from alternative splicing events. The splicing events included in the PI may provide promising therapeutic targets for HCC patients.

As mentioned above, alternative splicing events do not occur randomly, but they are rather regulated by multiple cellular factors, including splicing factors (56,57). Splicing factors can activate or inhibit the occurrence of splicing events by identifying regulatory sequences in exons or introns. DNA mutations, DNA methylation and aberrant histone modifications may all affect the recognition of splicing sites by splicing factors, leading to changes in the occurrence of splicing events (10,58-67). The roles of common anomalies, such as DNA mutations, in cancer can be rationally explained by selectivity. Mutations within introns and synonymous mutations in exons, which were previously believed not to affect the function of genes, have gradually been found to affect gene function by affecting splicing events (10,58,60). Splicing factors play a crucial regulatory role in the occurrence of splicing events, and mutations in their sequences or changes in their expression levels may affect splicing events (21). Changes in the expression levels of splicing factors are usually accompanied by changes in splicing events, but the pattern is not yet clear. Both the upregulation and downregulation of their expression can cause cancer or inhibit tumor growth by altering splicing events (25). Analysis of the association between splicing factors and survival in patients with HCC revealed that the expression levels of 32 splicing factors were correlated with survival time, with high expression being associated with shorter survival. It was hypothesized that the expression levels of these factors may be regulated by the similar events. For example, such an effect of the MYC oncogene has been reported (22). To elucidate the mechanism underlying the effect of splicing factors on the survival of patients with HCC, we analyzed the correlations between the expression levels of these 32 splicing factors and survival-related alternative splicing events. Surprisingly, the alternative splicing events that were positively correlated with survival were also positively correlated with the expression levels of multiple splicing factors, and most alternative splicing events that were negatively correlated with survival were also negatively correlated with the expression levels of the splicing factors. This finding indicates that the increase in the expression of multiple splicing factors can affect the survival of patients with HCC by synergistically regulating the alternative splicing events of genes.

However, there were several limitations to the present study. First, this study was based on an individual source of data from TCGA, without validation from independent cohorts, which will be carried out by our group in the future with clinical samples. Furthermore, the biological roles of the splicing events require further validation. In addition, the present study mainly focused on the splicing status of mRNAs. In fact, the splicing events of other types of RNA, including circular RNAs, also play important roles in the processes of tumors (68). Further research and effective arithmetic are required to broaden its application range (69). Finally, classic splicing factors were utilized in the present study, whereas it is possible that other RNA-binding proteins may also affect splicing events. Hence, the association between splicing events and a complete repertoire of RNA-binding proteins from established sources in the human genome should be further investigated.

Our study preliminarily demonstrated that alternative splicing events and the expression levels of splicing factors play important roles in the progression of HCC. On the basis of these results, it may be considered that splicing factors affect the occurrence and progression of HCC by regulating alternative splicing events. Some special alternative splicing events may affect the prognosis and progression of HCC by being regulated by their corresponding splicing factors. The constructed computational prognostic model based on HCC-specific alternative splicing events may be used as a molecular marker for the prognosis of HCC.

Currently, certain anticancer drugs targeting splicing factors, such as E7107 and FR901464 (29,70), have been developed, but they are not widely used in clinical practice due to their prominent side effects (21). The identified HCC-specific alternative splicing events may be used as molecular markers for the diagnosis and prognosis of HCC, and they can help physicians develop more precise targeted therapies. This approach has the potential to be widely applied in the field of HCC research.