The selection criteria for this study are as follows: (1) definite histological diagnosis of HCC and ICC; (2) definitive clinical data; (3) at least 30 days of overall survival after initial pathologic diagnosis (14–16); (4) complete RNA-sequencing data. Gene expression quantification data and related clinical data of the HCC and ICC were downloaded from TCGA database, and DESeq2 package was used to normalize the data portal (17, 18). Besides, the Percent Spliced In (PSI) value, which is a widely accepted indicator to quantify the AS events, was downloaded from TCGA SpliceSeq (19). To obtain the most reliable AS events set, we set a series of strict filter conditions (Percentage of samples with PSI value more than 0.75, average of PSI value more than 0.05) (20). UpSet plots were generated by the package of UpSetR (version 1.4.0) to display interactive sets between each types of AS events (21). In addition, Circos plots were generated by the software of Circos (version 0.69-6) to depict the details of splicing events and location of parent gene in whole chromosome (22).

To identify the tumor-specific splicing events between tumor tissue and normal tissue, the PSI value of patients was calculated (including 343 HCC tissues and 48 normal tissues, 33 ICC tissues and 8 normal tissues). Benjamini & Hochberg (BH) correction was used to adjust p-values. The AS events with the adj.p < 0.05 and|log2FC|>1 were considered to be significantly upregulated or downregulated. And the Venn diagram was developed to represent the differences between differentially expressed alternative splicing (DEAS) events and differentially expressed genetic (DEG). The parent genes of DEAS event were submitted to the String 11.0 online database for protein-protein interaction (PPI) analysis (23). The relationship network was then illustrated by Cytoscape (24). In addition, these parent genes were also used as the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses via Metascape.

Additionally, the “ConsensusClusterPlus” package was performed to classify patients based on the DEAS events (25). Subsequently, immune features were analyzed using ESTIMATE (26) and ssGSEA (27). The correlation analysis was conducted to clarify the relationship between DEAS clusters (or two cancer types) and immune characteristics (ESTIMATE Score, Stromal Score, Tumor Purity, Immune Score, Cytolytic activity, NK cells, etc.).

To standardize the PSI data, the median PSI value was used as a threshold to divide patients into two groups for each AS event. Univariate cox analysis and LASSO (alignment=lambda, nfold=10, gamma = c(0, 0.25, 0.5, 0.75, 1) analysis were used to identify potential prognostic factors (28). Later, the multivariate cox analysis was used to identify the independent risk factors for overall survival (OS) and disease-free survival (DFS) (named OS-DEAS events and DFS-DEAS events). And the risk scores of patients were calculated based on multivariate cox model (named OS-model and DFS-model, respectively).

Then, to combine the OS- or DFS-model with clinicopathological data, the nomograms were developed by rms package (v6.2.0) (29, 30). In addition, the area under the curve (AUC) of the receiver operating characteristic (ROC) and the consistency index (C-index) was calculated to evaluate the predictive ability of nomogram or other models (31).

The selection and cultivation of cell lines, the materials and methods of functional experiments, the collection and processing of tissue samples, and other methods can be found in the supplementary materials.

The data of splicing factors (SFs) which were validated in previous studies were downloaded from the SpliceAid 2 database (32). In addition, the expressions of SFs were downloaded from TCGA database, and the DESeq2 package was used to normalize the data portal (18). Correlation analyses were performed to determine the potential regulatory relationship between both OS- or DFS-DEAS events and SFs. In addition, crystal structures of SFs are obtained from Protein Databank.

All statistical analyses were conducted by R software (version 3.6.1) (24). Categorical data were performed using chi-square (χ2) test. Spearman’s rank correlation analysis was utilized for non-normal distribution data. Student’s t-test and ANOVA test were utilized to compare continuous variables. Survival curves were compared using log-rank test and performed using the Kaplan–Meier method. Pearson correlation was utilized for continuous variables that meet normal distribution. The results of Cox analysis were presented as the mean ± S.D., and P<0.05 (two-tailed) was considered statistically significant.

The pipeline of our research is shown in Figure 2 . A total of 343 HCC patients and 33 ICC patients were included in this study (the baseline characteristic of patients is listed in Supplementary Tables S1 , S2 ). Subsequently, 24,763 AS events in 8,434 genes were further detected in 343 HCC patients, and 28,147 AS events and 8,094 genes were detected in 33 ICC patients. These data indicated that one gene could have nearly three types of AS events. AS events include seven subtypes ( Supplementary Figure S1A ), which were all detected both in HCC and ICC ( Supplementary Figures S1B, C ), and ES was the most common AS type and ME is the rare type of AS in tumors. As shown in Supplementary Figures S1D, E , the UpSet plots of HCC and ICC showed the sets of each AS type. Moreover, two Circos-plots were developed to depict the details of AS events and location of parent gene in whole chromosome ( Supplementary Figures S1F, G ). The above results indicate that alternative splicing, which leads to the different arrangements and combinations of exons and introns, is responsible for the diversity of the transcriptome.

To identify the tumor-specific splicing events in tumor tissues and normal tissues, the comparison of PSI value between these tissues was performed. Finally, 384 DEAS events were identified from 336 genes in HCC. Meanwhile, in ICC, 749 DEAS events were found from 622 genes in ICC. The details of DEAS events are shown in Supplementary Tables S3 , S4 . Intriguingly, we found the AP type was the predominant DEAS mode in both HCC and ICC, and the significant differences in the distribution of seven splicing modes about DEAS events between HCC and ICC ( Supplementary Figure S2A ).

After the DEAS events of HCC and ICC were identified and the upregulation and downregulation DEAS events were displayed in the volcano plots ( Figures 3A, B ), unsupervised hierarchical consensus clustering was performed basing on DEAS events. The results showed that samples of cancer and normal tissues can be clearly separated into two groups, which means that the DEAS events identified above were convincing ( Figures 3C, D ). Moreover, we developed two Venn diagrams to depict the relationship of DEAS events and DEG in HCC ( Figure 3E ) or ICC ( Figure 3F ), and the details are shown in Supplementary Tables S5 , S6 . Intriguingly, whether in HCC or ICC, many genes (such as ADRA1A and KIF4A) displayed some opposite features of AS events in tumor and normal tissues ( Supplementary Figure S2B ). Moreover, 139 AS events and 144 DEG were identified as common DEAS events and DEG between HCC ( Figure 3G ) and ICC ( Figure 3H ).

Subsequently, the corresponding proteins of DEAS events were used to construct PPI networks for HCC ( Supplementary Figure S3A ) and ICC ( Supplementary Figure S3B ), and the PPI networks analyses demonstrated their interactive relationship in normal condition. We further filtered 10 hub genes from the protein-protein interaction network in HCC ( Supplementary Figure S3C ) and ICC ( Supplementary Figure S3D ), respectively. Two hub genes, Fibronectin 1 (FN1) and Serpin peptidase inhibitor clade A member 1 (SERPINA1), were identified as common genes between HCC and ICC, but other eight hub genes are different. Moreover, we analyzed the potential functions of DEAS events by GO and KEGG pathway ( Figures 3I, J ). The results suggested that GO categories related to the metabolic process, including “cofactor metabolic process” in HCC and “small molecule catabolic process” in ICC. And “cofactor metabolic process” and “lipid biosynthetic process” in GO analysis were identified as common metabolic processes between HCC and ICC. Moreover, KEGG pathways enriched were associated with tumorigenesis, such as “PI3K-Akt signaling pathway” in HCC and “Chemical carcinogenesis” in ICC. And only the “Ferroptosis” was confirmed as common KEGG pathway between HCC and ICC. Intriguingly, immune-related pathways were also enriched in HCC (not in ICC), such as “Complement and coagulation cascades” and “chemokine signaling pathway,” which indicated that DEAS events may be involved in immune microenvironment formation in HCC patients. These results prove that HCC and ICC are partly common in DEAS events and potential functional, and they may provide a reference for the study of CHCC. More importantly, more DEAS events (86%) and pathways (90%) are different between HCC and ICC.

Immune microenvironment is crucial to the development and recrudescence of tumors. Therefore, we explored the differences in immune checkpoints (33) ( Supplementary Table S7 ) and immune cell infiltration ( Supplementary Figure S4A ) between HCC and ICC and found that the immune microenvironment (65.9% of the immune checkpoints and 41.2% of the immune cells) between HCC and ICC is very different. More importantly, we are interested in the potential impact of alternative splicing on the immune microenvironment; thus, we performed a hierarchical consensus clustering analysis of HCC and ICC patients based on the hierarchical consensus clustering of DEAS events. Eventually, the HCC was divided into three groups ( Figure 4A and Supplementary Figures S4B – D ), and the ICC was divided into three groups ( Figure 4E and Supplementary Figures S4F – H ). We found that 45 immune checkpoints were significantly different between DEAS clusters in HCC, while only seven immune checkpoints in ICC ( Supplementary Tables S8 , S9 ). Interestingly, tryptophan 2,3-dioxygenase (TDO2) serves as a common immune metabolism checkpoint for HCC and ICC, which is consistent with the “cofactor metabolic process” previously identified as a common metabolic process between HCC and ICC. More importantly, PD-1 was significantly lower in cluster 2 than in cluster 1/3 ( Figure 4B ), while the expression of butyrophilin like 9 (BTNL9) in cluster 2 was significantly higher than that in cluster 1/3 in HCC ( Figure 4C ). In ICC, Poliovirus receptor (PVR), as another immunosuppression-related molecule, was significantly lower in cluster 2 than in clusters 1 and 3 ( Figure 4F ), while the expression of TNF receptor superfamily member 14 (TNFRSF14) in cluster 2 was significantly higher than that in cluster 3 ( Figure 4G ).

Next, we explored the differences in the infiltration of 34 types of immune cells in the DEAS clusters and found that 23 immune cells were considered significantly different between DEAS clusters in HCC ( Figure 4D ), while only one immune cell was considered significantly different between DEAS clusters in ICC ( Figure 4H ). In addition, the tumor microenvironment score and tumor stroma score based on DEAS clustering in HCC are meaningful, but not meaningful in ICC. Intriguingly, the components of “Th1 cells” show significant correlation with clusters both in HCC and ICC ( Supplementary Figures S4E, I ). In general, these results show that the immune microenvironment between HCC and ICC is very different, and more importantly, the correlation between AS and the immune microenvironment in HCC and ICC is extremely valuable for research.

To further evaluate the importance of AS for HCC and ICC, we used univariate survival analysis ( Supplementary Tables S10 – S11 ) and LASSO regression analysis ( Supplementary Figure S5 ) to analyze the impact of these DEAS on the survival and prognosis of HCC and ICC patients. And these DEAS events were selected in multivariate analysis. According to the results of multivariate survival analysis, five OS-DEAS events and six DFS-DEAS events were found as independent predictors of OS and DFS in patients with HCC ( Figure 5A ). In addition, in patients with ICC, three OS-DEAS events and three DFS-DEAS events were found as independent predictors ( Figure 5B ). Then we developed risk prediction formulas based on the above such DEAS events in HCC and ICC, respectively; the formulas are summarized in Supplementary Table S12 .

The ideal predictive model should consider the importance of clinical data. Univariate survival analysis ( Tables 1 , 2 ) and lasso regression ( Supplementary Figure S6 ) were performed to identify suitable clinical predictors. Afterwards, four nomograms based on risk scores of above OS or DFS models and clinical variables were developed to predict the OS and DFS in patients with HCC ( Figures 6A, C ) or ICC ( Figures 6E, G ). In addition, the corresponding calibration curves of nomograms showed good agreement between the probability of prediction and observation in 1-, 2-, and 3-years OS and DFS in HCC ( Figures 6B, D ) and in the 0.5-, 1-, and 2-year OS and DFS in ICC ( Figures 6F, H ). The c-index of nomogram was 0.732 (95% CI: 0.671–0.793) in OS-HCC group, 0.683 (95% CI: 0.624–0.742) in DFS-HCC group, 0.762 (95% CI: 0.649–0.875) in OS-ICC group, and 0.771 (95% CI: 0.621–0.921) in DFS-ICC group ( Table 3 ), and the ROC of four nomograms was also performed and is shown in Supplementary Figure S7 , respectively. We also validate nomograms internally by randomly drawing 70% of the original cohort, and the results show that the c-index of nomograms was 0.8 in OS-HCC group, 0.77 in DFS-HCC group, 0.86 in OS-ICC group, and 0.82 in DFS-ICC group. These results show that four nomograms have good stability and distinguishing ability. Meanwhile, the c-index and AUC of all single variables included in the four nomograms were also identified, and the results showed that the c-index and AUC of almost all single predictors were lower than the nomogram. Only when predicting the 1-year DFS of ICC patients is the predictive power of the nomogram model slightly lower than that of the DFS-model (Nomo AUC: 0.849 vs. DFS-model AUC: 0.853).

As in the previous study, we found that AA_ADHFE1_ID_084004 was a common independent risk event for both OS and DFS in HCC ( Figure 5A ; ADHFE1, Alcohol Dehydrogenase Iron Containing 1). In addition, AD_PIR_ID_008558 was also a common independent risk factor for both OS and DFS in ICC ( Figure 5B ; PIR, Pirin). For intuitively showing the differences of AA_ADHFE1_ID_084004 and AD_PIR_ID_008558 between tumor and normal tissues, we performed graphs in scatter plot ( Figures 5C, F ). Based on the median value of PSI, patients were divided in high-risk group and low-risk group. The Kaplan-Meier curves showed that there are significant differences in two groups ( Figures 5D, E, G, H ). It is proved again that AA_ADHFE1_ID_084004 and AD_PIR_ID_008558 may play an important role in the development and recrudescence of HCC and ICC, respectively.

To further verify that AA_ADHFE1_ID_084004 is very important to the development and recrudescence of HCC, we verify the function of two transcripts related to DEAS by experiments. The top left plot displays the splicing pattern of ADHFE1_203 and ADHFE1_207 [ADHFE1, through Alternate acceptor site (AA), an alternative splicing type, produced two transcripts] ( Figure 7A ), and the top right plot shows that the expression of ADHFE1_203 and ADHFE1_207 in HCC is significantly lower than that in adjacent non-tumor frozen tissues ( Figure 7B ). More importantly, MTT assay suggested that the overexpression of ADHFE1_203 and ADHFE1_207 inhibited the proliferation of HCC cells ( Supplementary Figure S8A and Figure 7C ), and Cell Cycle analysis suggested that the overexpression of ADHFE1_203 inhibited S phase, while the overexpression of ADHFE1_207 inhibited G1 phase ( Figure 7D ). In addition, we observed that the overexpression of ADHFE1_203 significantly inhibited migration and invasion in Huh7 cell line ( Figure 7E ). Intriguingly, MTT and transwell assays showed that compared with the overexpression of ADHFE1_207, the overexpression of ADHFE1_203 had a stronger ability to inhibit proliferation, invasion, and migration. These results directly indicated that DEAS events were important biological processes and had a potential clinical value.

So, what are the reasons for the emergence of important DEAS such as AA_ADHFE1_ID_084004 and AD_PIR_ID_008558? As we all know, SFs are important factors to regulate the DEAS events. Hence, we further studied which SF can regulate the production of OS- and DFS-DEAS events. Thus, correlation analyses between expression levels of the PSI values of these DEAS events and 71 SFs ( Supplementary Table S13 ) were conducted to explore the candidate regulation network in the HCC and ICC (34). As shown in Figures 8A, B , we can find that most of the SFs positively related with these DEAS events [67.11% (204/304) in HCC and 61.54% (16/26) in ICC]. In addition, we can also find that most single SF was correlated with more than one DEAS events, and the number of AS events correlated with some SFs even reach nine (DAZAP1 and HNRNPL in HCC). Further, we identified the SFs that significantly correlated with common DEAS events determined above (AA_ADHFE1_ID_084004 and AD_PIR_ID_008558). And we found that a lot of SFs correlated with AA_ADHFE1_ID_084004, but the T-cell intracellular antigen 1 (TIA1) is the most significantly correlated SFs (r=0.532, p<0.001). However, only SF Proline and Glutamine Rich (SFPQ) was identified correlated with AD_PIR_ID_008558 (r=−0.424, p=0.027).

More importantly, to better understand the specific splicing mechanism of SFs, we further explored the protein crystal structure of TIA1 and SFPQ. Analysis of TIA1 sequence and structures indicates it contains three RNA recognition motif (RRM), and each RRM structure consists of four antiparallel strands and two helices arranged in an alpha/beta sandwich, with beta sheet interacting with RNA molecules. Supplementary Figure S8B highlights interactions observed in the published crystal structures. Analysis of SFPQ sequences reveal it also contains two RRMs. Though available crystal structures of SFPQ do not contain RNA, structure-based alignment of SFPQ (PDB code 6NCQ) with TIA1 (PDB code 5O3J) indicates that SFPQ is capable of recognizing RNA ( Supplementary Figure S8C ).