Skip to main content

Immunological significance of prognostic alternative splicing signature in hepatocellular carcinoma

Abstract

Background

Hepatocellular carcinoma (HCC) ranks the sixth prevalent tumors with high mortality globally. Alternative splicing (AS) drives protein diversity, the imbalance of which might act an important factor in tumorigenesis. This study aimed to construct of AS-based prognostic signature and elucidate the role in tumor immune microenvironment (TIME) and immunotherapy in HCC.

Methods

Univariate Cox regression analysis was performed to determine the prognosis-related AS events and gene set enrichment analysis (GSEA) was employed for functional annotation, followed by the development of prognostic signatures using univariate Cox, LASSO and multivariate Cox regression. K-M survival analysis, proportional hazards model, and ROC curves were conducted to validate prognostic value. ESTIMATE R package, ssGSEA algorithm and CIBERSORT method and TIMER database exploration were performed to uncover the context of TIME in HCC. Quantitative real-time polymerase chain reaction was implemented to detect ZDHHC16 mRNA expression. Cytoscape software 3.8.0 were employed to visualize AS-splicing factors (SFs) regulatory networks.

Results

A total of 3294 AS events associated with survival of HCC patients were screened. Based on splicing subtypes, eight AS prognostic signature with robust prognostic predictive accuracy were constructed. Furthermore, quantitative prognostic nomogram was developed and exhibited robust validity in prognostic prediction. Besides, the consolidated signature was significantly correlated with TIME diversity and ICB-related genes. ZDHHC16 presented promising prospect as prognostic factor in HCC. Finally, the splicing regulatory network uncovered the potential functions of splicing factors (SFs).

Conclusion

Herein, exploration of AS patterns may provide novel and robust indicators (i.e., risk signature, prognostic nomogram, etc.,) for prognostic prediction of HCC. The AS-SF networks could open up new approach for investigation of potential regulatory mechanisms. And pivotal players of AS events in context of TIME and immunotherapy efficiency were revealed, contributing to clinical decision-making and personalized prognosis monitoring of HCC.

Introduction

Hepatocellular carcinoma (HCC) is a malignant and aggressive disease marked with frequently diagnosed and high cancer-attributable mortality in the world [1,2,3]. Despite great advance in HCC early diagnosis and anti-tumor therapy, clinical treatment result is still undesirable [4,5,6]. Both tumor, node, metastasis (TNM) categories and Barcelona Clinic Liver Cancer (BCLC) staging classification, which were widely used prognostic tools, failed to precisely predict results of patients with same clinicopathological stage [7,8,9]. Furthermore, the high heterogeneity of HCC remains great challenge against therapeutic benefits, which makes it difficult to accurately predict clinical result [1, 10, 11].

A great part of HCC was derived from inflammatory liver diseases, which suggested that infiltrating immune cells in tumor immune microenvironment (TIME) might serve as pivotal regulatory roles in tumorigenesis and progression in HCC [12,13,14]. Recently, immunotherapy has received extensive attention since it yielded encouraging results in multiple malignancies [15]. However, only 20% of HCC patients were observed objective response to immunotherapy according to preclinical trials [16]. Hence, the most effective strategy for precise prognostic predictions of how a given malignancy will respond to immunotherapy or tumor will progress may be one based on molecular risk distribution, identifying HCC patients on line with particular molecular signatures, enhancing prognostic precision and optimize immunotherapeutic benefit accordingly.

Alternative splicing (AS) defined as the mechanism by which edit pre-mRNA to produce mature mRNA, greatly contributed to the complexity of genome and the diversity of proteome [17, 18]. It is well established that AS events included such main patterns as alternate acceptor site (AA), alternate donor site (AD), alternate promoter (AP), alternate terminator (AT), exon skip (ES), retained intron (RI) and mutually exclusive exons (ME) [19]. Unbalanced expression of AS occurs frequently in cancer and received extensive attention as pivotal role in tumor initiation, development, metastasis and response to treatment [20,21,22]. Besides, splicing factors were found to act as a vital player in the regulation of AS events [23]. It was worth mentioned that aberrant expression of crucial splicing factors can lead to the oncogenic splicing isoforms [24, 25].To date, multiple researches had been performed to unravel the biological relevance and clinical significance of AS events in HCC [26,27,28].

There have been several articles focusing on AS-based prognostic model [29,30,31], however, the correlation of AS prognostic signature with TIME and immunotherapy remains obscure.

Thus, it is imperative to perform a comprehensive analysis of AS events to uncover the characterization of TIME and underlying molecular mechanisms of tumorigenesis, further optimize clinical benefits.

In this study, the AS pattern of TCGA LIHC cohort was outlined and AS events associated with the survival was determined through comprehensive bioinformatic analysis. Next, AS-based prognostic signatures were developed then validated, and an AS-clinicopathologic nomogram was constructed to facilitate clinical application. Then, the correlation of prognostic signature with the complexity of TIME and immunotherapeutic efficacy was investigated. Additionally, the potential role of ZDHHC16 in HCC was explored. Finally, an AS-SFs regulatory network was established to elucidate the potential mechanism involving in HCC progression.

Materials and methods

Multiomic data acquisition

The transcriptome information and survival information of the HCC patients were downloaded from The Cancer Genome Atlas (TCGA) portal (http://cancergenome.nih.gov) for subsequent analysis. The alternative splicing data of TCGA LIHC-cohort were obtained from SpliceSeq (http://bioinformatics.mdanderson.org/TCGASpliceSeq). Samples were screened when setting PSI value exceeds 0.75 as filter cut-off point. All analyses were performed based on the publication guidelines of TCGA. The analysis process flow chart was presented in Additonal file 1: Figure S1.

Process of AS profile identification

In TCGA splice-seq, the percent spliced in (PSI) value to quantify AS events were detected then calculated. By using UpsetR package, Upset plot was delineated to discovery the seven subtypes of AS events (alternate acceptor site (AA), alternate donor site (AD), alternate promoter (AP), alternate terminator (AT), exon skip (ES), mutually exclusive exons (ME), and retained intron (RI)). The AS events were annotated by combining the splicing type, ID number in the SpliceSeq and the corresponding parent gene symbol. For example, in “MRPL43|12,849|AT”, MRPL43 denotes the corresponding parent gene name, 12,849 represents the ID of splicing variant and AT indicates the splicing type.

Identification of survival-related AS events

When the standard deviation of PSI value is less than 0.01, the AS data were excluded. Univariate Cox regression analysis was carried out to detect the association between AS events and overall survival (OS) of HCC patients (Additonal file 2: Table S1), which were presented with the UpSet plot and the volcano map. Besides, the top 20 most significant AS events from the seven subtypes were summarized in the bubble charts.

Construction and validation of prognostic signature

Firstly, Lasso regression analysis was employed to determine candidates in each splicing pattern and avoid model over-fitting. Secondly, the identified AS events were introduced into Multivariate Cox regression analysis to screen the prognostic predictor. Given the mode of AS events is independent from each subtype in post-transcriptional modification, the identified AS events in each splicing subtype above were integrated to generate another prognostic signature. Subsequently, the risk scores were calculated according to each prognostic predictor and the formula for computing the risk score is as follows: Risk score = βAS event1 × PSIAS event1 + βAS event2 × PSIAS event2 +  + βAS eventn × PSIAS eventn. The specific formulas of each signature were presented in Additonal file 2: Table S2. Based on the median value of risk score, patients were ranked into low-risk group and high-risk group. Kaplan–Meier survival curves were analyzed with “survival” R package. Then, the time-dependent receiver operating characteristic (ROC) curves were performed to examine the prognostic value of this signature. Besides, univariate and multivariate Cox regression were analyzed to confirm whether the signature can serve as an independent factor for prognostic prediction. Additionally, stratified survival analysis was conducted to further validate the prognostic performance independent from such clinical characteristics as age, gender, tumor stage, pathological grade, T category, N category and M category.

Construction of prognostic nomogram

To comprehensively assess prognosis predictive ability of risk signature, tumor stage, gender, age, WHO grade, T category, N category and M category for 1/2/3-year OS, time-dependent receiver operating characteristic (ROC) curves was perform to calculate the area under the curve (AUC) values [32]. To contribute a quantitative manner to predicting overall survival of patients with HCC, prognostic nomogram which containing AS-based risk model and other clinical variables was established to estimate 1‐, 2‐and 3‐year overall survival probability. Subsequently, the calibration curve which shown the prognostic value of as-constructed nomogram was delineated. A calibration curve close to 45° is an indication of good prediction ability of the model constructed by this factor.

Correlation of risk score with tumor infiltrating immune cells characterization

Immune infiltration information consists of every specimen immune cell fraction (i.e., B cells, CD4 + T-cells, CD8 + T-cells, dendritic cells, macrophages, and neutrophils, etc.,) were downloaded from tumor immune estimation resource (TIMER) (https://cistrome.shinyapps.io/timer/). The correlation between tumor immune cell infiltrating with the prognostic risk score was performed. A single sample gene-set enrichment analysis (ssGSEA) was implemeted to elucidate the enrichment of the two distinct risky subgroups in 29 immune function‑associated gene sets via invoking the R package “GSEAbase”. Subsequently, the R package “ESTIMATE” was employed to assess tumor purity and the extent and level of infiltrating cells, namely stromal cell and immune cell, that could validate significant distinct tumor immune microenvironment characterization between two risky subgroups. The fraction of 22 immune cell types for each tumor specimen was developed through cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT; https://cibersort.stanford.edu/).

Role of risk score in immune checkpoint blockade treatment

Refer to existing studies, expression level of immune checkpoint blockade-related key genes might be correlated with clinical outcome of immune checkpoint inhibitors blockade treatment [33]. Herein, six key genes of immune checkpoint blockade therapy: programmed death ligand 1 (PD‐L1, also known as CD274), programmed death ligand 2 (PD‐L2, also known as PDCD1LG2), programmed death 1 (PD‐1, also known as PDCD1), cytotoxic T‐lymphocyte antigen 4 (CTLA‐4), indoleamine 2,3‐dioxygenase 1 (IDO1), and T‐cell immunoglobulin domain and mucin domain‐containing molecule‐3 (TIM‐3, also known as HAVCR2) in HCC [34,35,36] were extracted. To elucidate the potential player of as-constructed risk signature in ICB treatment of HCC, AS-based prognostic signature and expression level of six immune checkpoint blockade key genes were correlated. Finally, the expression level of 47 immune checkpoint blockade-related genes (i.e., PDCD1, etc.,) between low-/high-risk groups were compared.

Construction of splicing regulatory network

A list of 404 splicing factors (SFs, Additonal file 2: Table S3) was referred to a previous research [37] and the RNA-seq profiles of SFs were downloaded from the TCGA database. The Spearman correlation analysis was performed to evaluate the association between the SFs and the survival-relevant AS events (Additonal file 2: Table S4). P < 0.001 and Correlation coefficient > 0.6 was the cutoff values. Finally, Cytoscape (version 3.8.0) was employed to build an underlying SF-AS regulatory network.

Experimental validation

L02 cell (human hepatic cell line) and two human HCC cell lines (MHCC-97H cells and HCC-LM3) were purchased from the Cell Bank of the Type Culture Collection of the Chinese Academy of Sciences, Shanghai Institute of Biochemistry and Cell Biology. The cell lines were all cultured in Dulbecco’s minimum essential media (DMEM) plus 10% fetal bovine serum (FBS; Invitrogen, Carlsbad, CA, USA). All cell lines were grown without antibiotics in a humidified atmosphere of 5% CO2 and 99% relative humidity at 37℃. Three different cell lines were subjected to quantitative real-time polymerase chain reaction (qRT-PCR).

RNA isolation and qRT-PCR analysis

Total RNA was extracted from cells using TRIzol (Invitrogen, Carlsbad, CA, USA) according to provided instructions. RNA concentration and purity were measured in triplicates utilizing the NanoDrop 2000 spectrophotometer (Thermo Scientific Inc., Waltham, MA, 93 USA). Then, total RNA was reverse transcribed to cDNA using the cDNA Reverse Transcription Kit (Vazyme, Nanjing, China). To determine the expression of ZDHHC16, cDNAs were subjected to qRT-PCR using SYBR Green Real-time PCR Master Mix (Takara) in Applied Biosystems 7500/7500 Fast Real-Time PCR System (Thermo Fisher Scientific). All samples were analyzed in triplicates. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) levels were used as the endogenous control and relative expression of ZDHHC16 was calculated using the \({2}^{{ - \Delta \Delta {\text{C}}_{{\text{t}}} }}\) method. The sequences of primers used for PCR were as follows: ZDHHC16, 5′-CCACCAGACTCCACCACCTACC -3′ (forward) and 5′-GCCACAGAACTGCACAGGAACC -3′ (reverse); and GAPDH, 5′-CAGGAGGCATTGCTGATGAT-3′ (forward) and 5′-GAAGGCTGGGGCTCATTT-3′ (reverse).

Statistical analysis

The Wilcoxon test was employed to compare two groups, whereas the Kruskal–Wallis test was carried out to compare more than two groups. Overall survival (OS) refers to the interval from the date of diagnosis to the date of death. Survival curves were plotted via the Kaplan–Meier log rank test. Risk scores, clinical variables, immune cell infiltrating extent and immune checkpoints were correlated with Pearson correlation test. CIBERSORT algorithm results with p >  = 0.05 were rejected for further analysis. Univariate and multivariate analyses were performed via Cox regression models to validate the independent prognosis predictive performance of risk signature. The prognostic value of the AS-based signatures for 1-, 2- and 3-year OS was assessed with the ROC curves. p < 0.05 deemed as statistical significance. R software (version 4.0.2) was utilized for all statistical analyses.

Results

Clinical characteristics and integrated AS events profiles in HCC

377 HCC patients were obtained using the TCGA database, and seven patients with incomplete information were excluded from this study. In total, 370 patients were enrolled. The basic clinical information of patients is shown in Table 1. The AS events profiles were comprehensively analyzed and the gene intersections among the seven subtypes of AS events was presented in UpSet plot (Additonal file 1: Figure S2A). These results showed that ES was the predominant splicing pattern meanwhile AD marked as the least frequent.

Table 1 Baseline data of all HCC patients

Identification of the survival-relevant AS events

With the help of Univariate Cox regression analysis, a total of 3294 AS events which were significantly related with survival were identified (p < 0.05). The detailed description was recorded in TAdditonal file 2: able S1. The intersecting sets of genes and splicing subtypes were delineated in the UpSet plot (Additonal file 1: Figure S2B). Among these subtypes of AS events, ES was the predominant pattern. The volcano map was generated to display the AS events (Fig. 1a). The first 20 significant survival-related AS events from the seven subtypes were summarized in Figs. 1b–h, 2.

Fig. 1
figure1

The survival-relevant AS events. a) The volcano plots of survival-relevant AS events. The most significant survival-relevant AAs, ADs, APs, ATs, ESs, MEs and RIs in TCGA LIHC cohort (bh)

Fig. 2
figure2

Confirmation of ALL AS-based prognostic signature. a LASSO coefficient profiles of the whole AS events. b Ten‐time cross‐validation for tuning parameter selection in the lasso regression. c Heatmap of the ALL signature AS events PSI value in HCC. The color from red to blue shows a trend from high expression to low expression. d Distribution of ALL signature risk score. e The survival status and duration of HCC patients. f Kaplan–Meier curve presenting survival in the high-risk and low-risk sets. g ROC analysis of the risk scores for overall survival prediction. The AUC was calculated for ROC curves, and sensitivity and specificity were calculated to assess score performance. Proportional hazards model results. h Univariate Cox regression results. i Multivariate Cox regression results

Development of the prognostic signature

Stepwise Lasso algorithm and multivariate Cox regression analysis were employed to estimate the prognostic performance of these identified survival relevant AS events. The results of Lasso regression analysis including seven subtypes of AS events and amalgamated AS events with seven splicing types were displayed in Additonal file 1: Figures S3A–S3G, S4A–S4G and 2A–2B. Then, multivariate Cox analysis was implemented to determine optimal survival-relevant AS events. Lastly, eight AS (AA, AD, AP, AT, ES, ME, RI, and ALL) prognostic signature were constructed. Table 2 presented the detailed formulas of each signature.

Confirmation of the prognostic signature

Based on the cut-off value of median risk score, HCC patients were stratified into low and high-risk subgroups for further research. Heatmaps displayed the distributions of AS events PSI values with corresponding subgroups and patients (Additonal file 1: Figures S5A, S5D, S6A, S6D, S7A, S7D, S8A and 2C). The allocations of risk score (Additonal file 1: Figures S5B, S5E, S6B, S6E, S7B, S7E, S8B and 2D) and dot pot of survival status (Additonal file 1: Figures S5C, S5F, S6C, S6F, S7C, S7F, S8C and 2E) suggested that high-risk HCC patients had shorter overall survival. Besides, Kaplan–Meier curve corroborated that patients with low-risk possessed significant better prognosis than patients in high-risk group (Additonal file 1: Figures S9A, S9C, S9E, S9G, S10A, S10C, S10E and 2F; all P < 0.05). To assess the prognostic value of risk signatures in HCC cohort, ROC curve were further analyzed. Area under curves of risk scores at 1-, 2- and 3‐year survival times were all more than 0.75, suggesting great sensitivity and specificity of their survival predictive ability (Additonal file 1: Figures S9B, S9D, S9F, S9H, S10B, S10D, S10F and 2G). Besides, results of univariate Cox model (Additonal file 1: Figures S11A, S11B, S11C, S11D, S11I, S11J, S11K and 2H) and multivariate Cox regression analysis (Additonal file 1: Figures S11E, S11F, S11G, S11H, S11L, S11M, S11N and 2I), suggesting risk scores could act an independent indicator in HCC.

A stratification analysis was employed to validate whether ALL prognostic signature still had powerful prognostic predictive ability when HCC patients classified into various subgroups based on clinical characteristics. Relative to patients with low-risk, high-risk HCC patients presented poorer prognosis in both the early- and late-stage subgroups (Additonal file 1: Figures S12A, S12B). Similarly, ALL prognostic signature presented excellent prognostic prediction performance for patients in T1-2 or T3-4 status (Additonal file 1: Figures S12C and S12D), patients male or female gendered (Additonal file 1: Figures S12E, S12F), patients in 1–2 or 3–4 tumor grade (Additonal file 1: Figures S12G, S12H), patients aged <  = 65 years or > 65 years (FigAdditonal file 1: ures S12I, S12J), patients in N0 status (Additonal file 1: Figures S12K), and patients in M0 status (Additonal file 1: Figures S12L). These results suggested that it can be an outstanding predictor independent from clinical parameters in patients with HCC.

Correlation of ALL prognostic signature with clinical features and construction of AS-clinicopathological nomogram

Differences of risk score among different subtypes according to clinical variables were determined to uncover its clinical significance. The risk score increased significantly with advanced tumor grade (most p < 0.05, Fig. 3a), advanced clinicopathological stage (most p < 0.05, Fig. 3b) and advanced T category (most p < 0.05, Fig. 3c), suggesting risk score was positively correlated with tumor progression. To explore whether ALL prognostic signature was the best prognostic indicator among various clinical characteristics, age, gender, clinicopathological stage, and tumor grade were extracted as the candidate prognosis predictive factors. These clinical variables were consolidated to conduct the AUC curve analysis for 1-, 2-, and 3-year OS and risk signature obtained the most AUC value (Fig. 3d–f). Then, prognostic nomogram including risk score and clinicopathological stage were established to forecast prognosis of patients with HCC (Fig. 3g). Age, gender and tumor grade were rejected out of the nomogram because of their AUCs were less than 0.6. Calibrate curves were approximately diagonal, indicating powerful prognostic predictive ability of 1-, 2- and 3-year OS in our nomogram plot (Fig. 3h–j).

Fig. 3
figure3

Correlation of ALL prognostic signature with clinical features and construction of AS-clinicopathological nomogram. a Correlation of risk score with tumor grade. b Correlation of risk score with clinicopathological. c Correlation of risk score with T status. df Areas under curves (AUCs) for predicting 1-, 2-, and 3-year survival with different clinical characteristics. g Nomogram was assembled by stage and risk signature for predicting survival of HCC patients. h One‐year nomogram calibration curves. i Two‐year nomogram calibration curves. j Three‐year nomogram calibration curves

Correlation of risk score with tumor immune environment characterization

To further examine whether risk score can act as immune indicators, the correlation analysis of prognostic risk score with TICs from TIMER, immune score (calculated using the ESTIMATE algorithm), ssGSEA signatures and TICs subtype and level (calculated via CIBERSORT method) were carried out. Firstly, TIMER results showed that the as-constructed signature exhibited the marked positive association with B cells infiltration (r = 0.116; p = 0.026), CD8 + T cells infiltration (r = 0.223; p = 2.349e − 05), Dendritic cells infiltration (r = 0.228; p = 1.564e − 05), Macrophages infiltration (r = 0.271; p = 2.357e − 07) and Neutrophils infiltration (r = 0.221; p = 2.945e − 05; Fig. 4a–e), indicating high-risk samples were infiltrated more immune cells. Likewise, low-risk patients obtained a higher stromal score which represented less immune infiltration (Fig. 4f). However, there was no significant difference regarding to immune score, estimate score and tumor purity (Additonal file 1: Figures S13A, S13B and S13C). Subsequently, distinction of the immune-related signatures between these two subgroups was presented. Figure 4g and h showed that immune-related signature of each patient with corresponding immune scores in low-/high-risk groups. The results showed that the infiltration of aDCs, Macrophages, Neutrophils, Tfh, Th1 cells, Th2 cells, Tregs, HLA molecule expression level and MCH class I expression, such immune signatures as APC costimulation, T cells costimulation, check-point, inflammation-promoting, and IFN response were significantly increased with increased risk score (Fig. 4i). The CIBERSORT algorithm results indicated that proportion of reseing Dendritic cells was negatively associated with risk score (Fig. 4j). Above results indicated that ALL prognostic signature may provide a novel approach to elucidate the characteristics of immunity regulatory network in HCC.

Fig. 4
figure4

Correlation between infiltrating immune cells and ALL AS-based prognostic signature. a Relationship between this signature and B cells. b Relationship between this signature and CD8 + T cells. c Relationship between this signature and Dendritic cells. d Relationship between this signature and Macrophages. e Relationship between this signature and Neutrophils. (F) Comparison of stromal score between low- and high-risk groups. g Heatmap displayed enrichment of 29 immune signatures of low-/high-risk groups. Blue represents low activity and red represent high activity. h Heatmap of 29 immune signatures and immune scores of two different risk score groups. Blue represents low activity and red represent high activity. i Distinction of enrichment of immune-related signatures between risk-low and risk-high groups. j Difference of infiltrating immune cell subpopulations and levels between low-/high-risk groups

Correlation of ALL signature with ICB key molecules

With the emergence of immune checkpoint blockade (ICB) therapy, immune checkpoint inhibitors have considerably transformed clinical decision-making in cancer oncology [38,39,40]. Subsequently, six key immune checkpoint inhibitors genes (PDCD1, CD274, PDCD1LG2, CTLA‐4, HAVCR2, and IDO1) [34,35,36] were correlated. And the correlation between ICB key targets and ALL prognostic signature was analyzed to reveal the potential player of risk signature in the ICB treatment of HCC (Fig. 5a). The results indicated that ALL prognostic signature was significantly positive correlated to CD274 (r = 0.26; p = 0.00015), CTLA4 (r = 0.33; p = 1.3e − 06), HAVCR2 (r = 0.41; p = 1.4e − 09), IDO1 (r = 0.15; p = 0.03), PDCD1 (r = 0.16; p = 0.021) and PDCD1LG2 (r = 0.23; p = 0.001; Fig. 5b–g). Further correlation analysis presented that 33 of 47 (i.e., PDCD1, CTLA4, etc.,) immune check blockade-associated genes expression levels were significantly upregulated in patients with high-risk (Fig. 5h), suggesting ALL prognostic signature might act as nonnegligible and unfavorable factor in immunotherapy operating.

Fig. 5
figure5

Association between ALL AS-based prognostic signature and crucial immune checkpoint genes. a association analyses between immune checkpoint inhibitors CD274, PDCD1, PDCD1LG2, CTLA4, HAVCR2, and IDO1 and risk score. b association between risk score and CD274. c association between risk score and CTLA4. d association between risk score and HAVCR2. e association between risk score and IDO1. f association between risk score and PDCD1. g association between risk score and PDCD1LG2. h Comparison of immune checkpoint blockade-related genes expression levels between low-risk group and high-risk groups

ZDHHC16 independently affected prognosis and correlated with ICB key genes

ZDHHC16 was only one gene whose expression level was upregulated among the prognostic AS-related genes. Therefore, the role of ZDHHC16 in HCC was explored in further experimental validation. ZDHHC16 expression level between normal tissues and tumor samples was compared based on TCGA data. Relative to tumor tissues, ZDHHC16 expression level was lower in adjacent normal specimens (Fig. 6a). Taking advantage of qRT-PCR, ZDHHC16 expression level in two distinct HCC cell lines and human hepatic cell line were determined. Consistent of results of online database, ZDHHC16 was upregulated in cancer cells relative to normal cell (Fig. 6b). The expression level analysis among major pathological stages presented that ZDHHC16 expressed statistical significantly among different pathological stages (Fig. 6c, F = 3.45 and P = 0.0168). Additionally, the later tumor grade, the higher risk score (Fig. 6d, almost p < 0.05). To further assess the prognostic value of ZDHHC16 in HCC, Kaplan–Meier analysis were conducted between ZDHHC16 low- and high-expressed patients. As presented in Fig. 6e and f, lower ZDHHC16 expression level significantly suggested longer overall survival time (p = 0.0056) and longer disease-free survival time (p = 0.02). Besides, 16 of 47 immune check blockade-associated genes (i.e., PDCD1, CTLA4, etc.,) expression levels between low-ZDHHC16 and high-ZDHHC16 groups were significantly dysregulated in between different subgroups (Fig. 6g). Then the correlation between the ZDHHC16 and ICB key targets adjusted by tumor purity using TIMER was analyzed to investigate the potential player of ZDHHC16 in ICB treatment of HCC. TIMER results presented ZDHHC16 was significantly positive correlated to CD274 (r = 0.132; p = 1.41e − 02), CTLA4 (r = 0.254; p = 1.79e − 06), HAVCR2 (r = 0.231; p = 1.50e − 05) and PDCD1 (r = 0.291; p = 3.66e − 08; Fig. 6h–k), suggesting ZDHHC16 may exert a vital player in ICB treatment of HCC.

Fig. 6
figure6

The clinical significance of ZDHHC16 in HCC and in vitro study. ZDHHC16 are overexpressed in HCC tumor tissue (a) and HCC cell lines (b). c The expression of ZDHHC16 had significant difference between major pathological stages. d Correlation of risk score with tumor grade. Lower ZDHHC16 level predicts longer overall survival (e) and disease-free survival (f). g Comparison of immune checkpoint blockade-related genes expression levels between low-ZDHHC16 group and high-ZDHHC16 group. h Correlation of risk score with CD274. i Correlation of risk score with CTLA4. j Correlation of risk score with HAVCR2. (K) Correlation of risk score with PDCD1

Role of ZDHHC16 in context of TIME

To further elucidate the relationship between ZDHHC16 and TIME characteristics in HCC, comprehensive analysis were performed as descripted previously. HCC patients were classified into high-/low- ZDHHC16 subtypes based on the median ZDHHC16 expression level. ESTIMATE results indicated that patients with lower ZDHHC16 expression had a significant higher stromal score and higher immune score relative to patients in high- ZDHHC16 group, suggesting less stromal cells and immune cells in low-risk samples. (Figs. 7a, b). Additionally, arm-level deletion was predominant type of mutation (Fig. 7c). Subsequently, expression level of ZDHHC16 was positively correlated with infiltration of main immune cells types (Fig. 7d). outcomes of ssGSEA showed that the infiltration fraction of B cells, neutrophils, NK cells, T helper cells and TIL, APC co-inhibition, T cell co-inhibition, CCR, cytolytic activity, IFN-response type-I and HLA expression were significantly increased when risk score declining (Fig. 7e). The CIBERSORT analysis results of TCGA cohort showed that the proportion of activated memory CD4 T cells was significantly higher in patients with low-risk (Fig. 7f).

Fig. 7
figure7

The role of ZDHHC16 in TIME features. a Comparison of stromal score, immune score and ESTIMATE score between low-/high-ZDHHC16 groups. b Comparison of tumor purity between low-/high-ZDHHC16 groups. c Copy number of immune cells in HCC. d Relationship between risk score with B cells, CD8 T cells, CD4 T cells, Macrophages, Neutrophils and Dendritic cells. e Comparison of ssGSEA enrichment between low-/high-ZDHHC16 groups. f Comparison of CIBERSORT results between low-/high-ZDHHC16 groups

Development of the SF-AS regulatory network

To elucidate the underlying mechanism of AS regulation, a correlation network between the expression level of SFs and the PSI values of prognosis-related AS events was constructed. 55 up-regulated AS events (yellow ellipses), 56 down-regulated AS events (blue ellipses) and 106 SFs (Fig. 8) were identified. In our regulatory network, the top four most significant nodes were termed as the hub SFs or AS events (Additonal file 2: Table S4), including one downregulated AS event (ACAA1-64,022-ES), one upregulated AS events (SCP2-3045-ES) and two SFs (ISY1 and CLK2). As such, these SFs exhibited promising potential to serve as pivotal regulators involving in the dysregulation of AS in HCC, further mediated tumor initiation and progression.

Fig. 8
figure8

The regulatory network between SFs and survival related AS events. The yellow or blue ellipses indicated the AS events positively or negatively correlated with survival. Purple rectangles represented SFs. The positive/negative correlations (r > 0.6 or r < − 0.6) between SFs and AS events were indicated with red/green lines

Discussion

As one of the most common type of malignant tumor, hepatocellular carcinoma (HCC) had high cancer-relevant mortality globally [1,2,3]. Since such intricate molecular mechanism as genomic complexities and epigenetic diversities, HCC was highly heterogeneous from both clinical standpoint and molecular level [41,42,43]. A sober reality is that a majority of HCC patients cannot obtain benefit from immunotherapy, due to tumor-promoting condition mediated by immunosuppressive cells (i.e., regulatory T cells, etc.,) [44]. Thus, there is an urgent call to develop powerful prognostic tools for immunotherapeutic outcome prediction, which could contribute novel insight into individual tailored treatment in HCC.

Increasing studies have provided strong evidence to support that AS, which refers to post‑transcriptional modification procedure, function in physiological and pathological process [17]. Notably, abnormally regulated AS events participated well in tumor initiation and development, including HCC [21, 45]. Furthermore, dysregulated expressed genes can be employed as novel prognostic indicator and promising therapeutic targets. However, little to know the correlation of AS events with context of TIME and immunotherapeutic results in HCC.

In current study, AS data was obtained from TCGA SpliceSeq and comprehensive analysis of AS events in HCC samples. Taking advantage of univariate Cox regression analysis, 3294 AS events significantly associated with the survival were identified to explore the prognostic value of AS events. Next, prognostic signatures for HCC patients were proposed using systematic bioinformatic analysis. All eight (AA, AD, AP, AT, ES, ME, RI, and ALL) prognostic predictive signatures constructed by AS patterns presented powerful capability for the prognostic prediction in HCC. Notably, these AS-based prognostic signatures were robustly demonstrated by K-M survival analysis, ROC curve and Cox regression analysis. Furthermore, this signature retained excellent prognostic performance when HCC cases divided into groups based on clinicopathological factors. To transform ALL risk model into further clinical practice, a nomogram graph including prognostic signature with clinicopathological stage was plotted, and there was high consistence between predicted outcome and actual outcome. Besides, the most significant associated SF-AS regulatory network in the TCGA LIHC cohort were screened.

To reveal the role of AS events in the context of TIME in HCC, TIMER database, ESTIMATE algorithm, ssGSEA method and CIBERSORT analysis were conducted. These results presented that high-risk score group was marked with high infiltration of immune cells and more activated immune condition, which might promote immune recognition and trigger anti-tumor effect. These outcomes implied that risk score might facilitate immunotherapy results prediction. Since no ICB treatment dataset in HCC cohort, it was unable to investigate the relationship between risk score and ICB immunotherapy response. Then, risk score was positively and significantly correlated with six ICB key targets (i.e., CD274 and CTLA4) and 33 (i.e., PDCD1LG2, etc.,) immune check blockade-associated genes expression levels, which imply that risk score might contribute to strategizing the tailored immunotherapy.

ZDHHC16 was a DHHC encoding protein, which was tightly correlated with protein palmitoylation in previous researches [46]. Wei Shi et al. reported that ZDHHC16 acted as a vital regulator in the process of NSPCs proliferation [47]. Li Jian et al. uncovered pivotal players for ZDHHC16 in the regulation of DNA damage responses and Atm activation(48). To date, little to know about the role of ZDHHC16 in tumors, especially in HCC. This study presents that ZDHHC16 is significantly upregulated in HCC cell lines and suggested poor prognosis in HCC. ZDHHC16 expression was significantly positive associated with clinicopathological stage, tumor grade and ICB immunotherapy key genes (i.e., CD274, CTLA4, HAVCR2 and PDCD1, etc.). Nevertheless, further research needed to explore the underlying biological roles of ZDHHC16.

Collectively, subjects with higher risk score or higher ZDHHC16 expression level presented more abundance of immune cells in tumor environment, suggesting activated immune phenotype, but shorter overall survival time. It could be speculated that anti-tumor effect of immune cells may be affected by immune checkpoint blockade pathways given risk score was correlated with immune checkpoint blockade targets expression.

Compared this research with existed studies that explored the novel prognostic factor in HCC, some superiorities of our research should be noted. Firstly, this study contributed to investigate the potential roles of AS events in formation of TIME diversity and complexity and ICB treatment prediction, which has not been elucidated before this study. Additionally, ESTIMATE R package, ssGSEA algorithm and CIBERSORT method and TIMER database exploration were performed to uncover the comprehensive landscape of TIME in HCC. Finally, to our knowledge, this work is the first placed emphasis on the biological functions of ZDHHC16 in HCC.

Conclusion

Collectively, systematical analyses in prognosis predictive value of RNA splicing patterns were performed, which was designed to strengthen prognosis prediction in HCC. It is worthwhile mentioned that novel and robust prognostic nomogram to predict outcome quantitatively was established, which exhibited encouraging potential into clinical application. Besides, the AS-SFs regulatory network suggested promising targets of the anti-tumor therapy in HCC. The comprehensive bioinformatic analysis of AS events robustly linked the AS atlas with TIME characterization and immunotherapy in HCC. Nevertheless, these findings should be validated in further experimental and clinical exploration which focusing on HCC tumorigenesis and progression mechanisms and the impacts of these AS events.

Availability of data and materials

The study was based on the data available at TCGA(https://www.cancer.gov/tcga).

Abbreviations

AA:

Alternate acceptor site

AD:

Alternate donor site

AP:

Alternate promoter

AT:

Alternate terminator

AUC:

Area under the curve

AS:

Alternative splicing

BCLC:

Barcelona Clinic Liver Cancer

CTLA‐4:

Cytotoxic T‐lymphocyte antigen 4

CD274:

Also known as PD-L1

DMEM:

Dulbecco’s minimum essential media

ES:

Exon skip

FBS:

Fetal bovine serum

GAPDH:

Glyceraldehyde-3-phosphate dehydrogenase

GO:

Gene ontology

HCC:

Hepatocellular carcinoma

HAVCR2:

Also known as TIM3

IDO1:

Indoleamine 2,3‐dioxygenase 1

ICB:

Immune checkpoint blockade

LASSO:

Least absolute shrinkage and selection operator

ME:

Mutually exclusive exons

OS:

Overall survival

PD‐1:

Programmed Cell Death 1

PD‐L1:

Programmed Cell Death-Ligand 1

PD‐L2:

Programmed Cell Death-Ligand 2

PDCD1:

Also known as PD-1

PDCD1LG2:

Also known as PD‐L2

qRT-PCR:

Quantitative real-time polymerase chain reaction

RI:

Retained intron

RNA:

Ribonucleic Acid

ROC:

Receiver operating characteristic

SFs:

Splicing factors

ssGSEA:

Single-sample gene set enrichment analysis

TCGA:

The Cancer Genome Atlas

TICs:

Tumor‐infiltrating immune cells

TIME:

Tumor immune microenvironment

TIMER:

Tumor immune estimation resource

TIM‐3:

T‐cell immunoglobulin domain and mucin domain‐containing molecule‐3

TNM:

Tumor, node, metastasis

References

  1. 1.

    Forner A, Reig M, Bruix J. Hepatocellular carcinoma. Lancet. 2018;391(10127):1301–14. https://doi.org/10.1016/s0140-6736(18)30010-2.

    Article  PubMed  Google Scholar 

  2. 2.

    Bray F, Ferlay J, Soerjomataram I, Siegel R, Torre L, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2018;68(6):394–424. Doi: https://doi.org/10.3322/caac.21492.

  3. 3.

    Yang J, Hainaut P, Gores G, Amadou A, Plymoth A, Roberts L. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat Rev Gastroenterol Hepatol. 2019;16(10):589–604. https://doi.org/10.1038/s41575-019-0186-y.

    Article  PubMed  PubMed Central  Google Scholar 

  4. 4.

    Song P, Cai Y, Tang H, Li C, Huang J. The clinical management of hepatocellular carcinoma worldwide: A concise review and comparison of current guidelines from 2001 to 2017. Biosci Trends. 2017;11(4):389–98. https://doi.org/10.5582/bst.2017.01202.

    Article  PubMed  Google Scholar 

  5. 5.

    Grandhi MS, Kim AK, Ronnekleiv-Kelly SM, Kamel IR, Ghasebeh MA, Pawlik TM. Hepatocellular carcinoma: From diagnosis to treatment. Surg Oncol. 2016;25(2):74–85. https://doi.org/10.1016/j.suronc.2016.03.002.

    Article  PubMed  Google Scholar 

  6. 6.

    Hambardzumyan M, Hayrapetyan A. Differential diagnosis of malignant melanoma and benign cutaneous lesions by ultrasound analysis. SciMed J. 2020;2:100–7. https://doi.org/10.28991/SciMedJ-2020-0202-7.

    Article  Google Scholar 

  7. 7.

    Liu P, Hsu C, Hsia C, Lee Y, Su C, Huang Y, et al. Prognosis of hepatocellular carcinoma: Assessment of eleven staging systems. J Hepatol. 2016;64(3):601–8. https://doi.org/10.1016/j.jhep.2015.10.029.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Chen L, Chang Y, Chang Y. Survival Predictability Between the American Joint Committee on Cancer 8th Edition Staging System and the Barcelona Clinic Liver Cancer Classification in Patients with Hepatocellular Carcinoma. Oncologist 2020 Doi: https://doi.org/10.1002/onco.13535.

  9. 9.

    Kosvyra A, Maramis C, Chouvarda I. Developing an integrated genomic profile for cancer patients with the use of NGS data. Emerg Sci J. 2019;3:157–67. https://doi.org/10.28991/esj-2019-01178.

    Article  Google Scholar 

  10. 10.

    Zhang Q, Lou Y, Yang J, Wang J, Feng J, Zhao Y, et al. Integrated multiomic analysis reveals comprehensive tumour heterogeneity and novel immunophenotypic classification in hepatocellular carcinomas. Gut. 2019;68(11):2019–31. https://doi.org/10.1136/gutjnl-2019-318912.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  11. 11.

    Labh S. Effects of lapsi choerospondias axillaris on growth and immune-related genes in silver carp (Hypophthalmichthys Molitrix). SciMed J. 2020;2:86–99. https://doi.org/10.28991/SciMedJ-2020-0202-6.

    CAS  Article  Google Scholar 

  12. 12.

    Ringelhan M, Pfister D, O’Connor T, Pikarsky E, Heikenwalder M. The immunology of hepatocellular carcinoma. Nat Immunol. 2018;19(3):222–32. https://doi.org/10.1038/s41590-018-0044-z.

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Nishida N, Kudo M. Immunological microenvironment of hepatocellular carcinoma and its clinical implication. Oncology. 2017;92(1):40–9. https://doi.org/10.1159/000451015.

    Article  PubMed  Google Scholar 

  14. 14.

    Zare H. Effects of salvia officinalis extract on the breast cancer cell line. SciMed J. 2019;1:25–9. https://doi.org/10.28991/SciMedJ-2019-0101-4.

    Article  Google Scholar 

  15. 15.

    Grosser R, Cherkassky L, Chintala N, Adusumilli PS. Combination immunotherapy with CAR T cells and checkpoint blockade for the treatment of solid tumors. Cancer Cell. 2019;36(5):471–82. https://doi.org/10.1016/j.ccell.2019.09.006.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  16. 16.

    Cheng H, Sun G, Chen H, Li Y, Han Z, Li Y, et al. Trends in the treatment of advanced hepatocellular carcinoma: immune checkpoint blockade immunotherapy and related combination therapies. Am J Cancer Res. 2019;9(8):1536–45.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. 17.

    Nilsen T, Graveley B. Expansion of the eukaryotic proteome by alternative splicing. Nature. 2010;463(7280):457–63. https://doi.org/10.1038/nature08909.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  18. 18.

    Black D. Mechanisms of alternative pre-messenger RNA splicing. Annu Rev Biochem. 2003;72:291–336. https://doi.org/10.1146/annurev.biochem.72.121801.161720.

    CAS  Article  PubMed  Google Scholar 

  19. 19.

    Li Y, Sun N, Lu Z, Sun S, Huang J, Chen Z, et al. Prognostic alternative mRNA splicing signature in non-small cell lung cancer. Cancer Lett. 2017;393:40–51. https://doi.org/10.1016/j.canlet.2017.02.016.

    CAS  Article  PubMed  Google Scholar 

  20. 20.

    Calabrese C, Davidson N, Demircioğlu D, Fonseca N, He Y, Kahles A, et al. Genomic basis for RNA alterations in cancer. Nature. 2020;578(7793):129–36. https://doi.org/10.1038/s41586-020-1970-0.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    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. https://doi.org/10.1016/j.celrep.2017.08.012.

    CAS  Article  PubMed  Google Scholar 

  22. 22.

    Lee SC, Abdel-Wahab O. Therapeutic targeting of splicing in cancer. Nat Med. 2016;22(9):976–86. https://doi.org/10.1038/nm.4165.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Tripathi V, Ellis J, Shen Z, Song D, Pan Q, Watt A, et al. The nuclear-retained noncoding RNA MALAT1 regulates alternative splicing by modulating SR splicing factor phosphorylation. Mol Cell. 2010;39(6):925–38. https://doi.org/10.1016/j.molcel.2010.08.011.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  24. 24.

    El Marabti E, Younis I. The cancer spliceome: reprograming of alternative splicing in cancer. Front Mol Biosci. 2018;5:80. https://doi.org/10.3389/fmolb.2018.00080.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  25. 25.

    Yang Q, Zhao J, Zhang W, Chen D, Wang Y. Aberrant alternative splicing in breast cancer. J Mol Cell Biol. 2019;11(10):920–9. https://doi.org/10.1093/jmcb/mjz033.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Li S, Hu Z, Zhao Y, Huang S, He X. Transcriptome-wide analysis reveals the landscape of aberrant alternative splicing events in liver cancer. Hepatology (Baltimore, MD). 2019;69(1):359–75. https://doi.org/10.1002/hep.30158.

    CAS  Article  Google Scholar 

  27. 27.

    Lee SE, Alcedo KP, Kim HJ, Snider NT. Alternative splicing in hepatocellular carcinoma. Cell Mol Gastroenterol Hepatol. 2020;10(4):699–712. https://doi.org/10.1016/j.jcmgh.2020.04.018.

    Article  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Zhu G, Zhou Y, Qiu L, Wang B, Yang Y, Liao W, et al. Prognostic alternative mRNA splicing signature in hepatocellular carcinoma: a study based on large-scale sequencing data. Carcinogenesis. 2019. https://doi.org/10.1093/carcin/bgz073.

    Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Cai Y, Xia J, Wang N, Zhou H. Identification of prognostic alternative splicing signatures in hepatitis B or/and C viruses related hepatocellular carcinoma. Genomics. 2020;112(5):3396–406. https://doi.org/10.1016/j.ygeno.2020.06.002.

    CAS  Article  PubMed  Google Scholar 

  30. 30.

    Wu F, Chen Q, Liu C, Duan X, Hu J, Liu J, et al. Profiles of prognostic alternative splicing signature in hepatocellular carcinoma. Cancer Med. 2020;9(6):2171–80. https://doi.org/10.1002/cam4.2875.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Zhang D, Duan Y, Wang Z, Lin J. Systematic profiling of a novel prognostic alternative splicing signature in hepatocellular carcinoma. Oncol Rep. 2019;42(6):2450–72. https://doi.org/10.3892/or.2019.7342.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  32. 32.

    Blanche P, Dartigues J, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med. 2013;32(30):5381–97. https://doi.org/10.1002/sim.5958.

    Article  PubMed  Google Scholar 

  33. 33.

    Goodman A, Patel S, Kurzrock R. PD-1-PD-L1 immune-checkpoint blockade in B-cell lymphomas. Nat Rev Clin Oncol. 2017;14(4):203–20. https://doi.org/10.1038/nrclinonc.2016.168.

    CAS  Article  PubMed  Google Scholar 

  34. 34.

    Kim J, Patel M, Mangraviti A, Kim E, Theodros D, Velarde E, et al. Combination therapy with Anti-PD-1, Anti-TIM-3, and focal radiation results in regression of murine gliomas. Clin Cancer Res. 2017;23(1):124–36. https://doi.org/10.1158/1078-0432.Ccr-15-1535.

    CAS  Article  PubMed  Google Scholar 

  35. 35.

    Nishino M, Ramaiya N, Hatabu H, Hodi F. Monitoring immune-checkpoint blockade: response evaluation and biomarker development. Nat Rev Clin Oncol. 2017;14(11):655–68. https://doi.org/10.1038/nrclinonc.2017.88.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Zhai L, Ladomersky E, Lenzen A, Nguyen B, Patel R, Lauing K, et al. IDO1 in cancer: a Gemini of immune checkpoints. Cell Mol Immunol. 2018;15(5):447–57. https://doi.org/10.1038/cmi.2017.143.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  37. 37.

    Seiler M, Peng S, Agrawal A, Palacino J, Teng T, Zhu P, et al. Somatic mutational landscape of splicing factor genes and their functional consequences across 33 cancer types. Cell Rep. 2018;23(1):282-96.e4. https://doi.org/10.1016/j.celrep.2018.01.088.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  38. 38.

    Llovet J, Montal R, Sia D, Finn R. Molecular therapies and precision medicine for hepatocellular carcinoma. Nat Rev Clin Oncol. 2018;15(10):599–616. https://doi.org/10.1038/s41571-018-0073-4.

    Article  PubMed  Google Scholar 

  39. 39.

    Pitt J, Vétizou M, Daillère R, Roberti M, Yamazaki T, Routy B, et al. Resistance mechanisms to immune-checkpoint blockade in cancer: tumor-intrinsic and -extrinsic factors. Immunity. 2016;44(6):1255–69. https://doi.org/10.1016/j.immuni.2016.06.001.

    CAS  Article  PubMed  Google Scholar 

  40. 40.

    Salik B, Smyth M, Nakamura K. Targeting immune checkpoints in hematological malignancies. J Hematol Oncol. 2020;13(1):111. https://doi.org/10.1186/s13045-020-00947-6.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Ally A, et al. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell. 2017;169(7):1327-41.e23. https://doi.org/10.1016/j.cell.2017.05.046.

    CAS  Article  Google Scholar 

  42. 42.

    Schulze K, Nault J, Villanueva A. Genetic profiling of hepatocellular carcinoma using next-generation sequencing. J Hepatol. 2016;65(5):1031–42. https://doi.org/10.1016/j.jhep.2016.05.035.

    CAS  Article  PubMed  Google Scholar 

  43. 43.

    Woo H, Kim Y. Multiplatform genomic roadmap of hepatocellular carcinoma: a matter of molecular heterogeneity. Hepatology (Baltimore, MD). 2018;68(5):2029–32. https://doi.org/10.1002/hep.29925.

    Article  Google Scholar 

  44. 44.

    Prieto J, Melero I, Sangro B. Immunological landscape and immunotherapy of hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. 2015;12(12):681–700. https://doi.org/10.1038/nrgastro.2015.173.

    CAS  Article  PubMed  Google Scholar 

  45. 45.

    Lu Y, Xu W, Ji J, Feng D, Sourbier C, Yang Y, et al. Alternative splicing of the cell fate determinant Numb in hepatocellular carcinoma. Hepatology (Baltimore, MD). 2015;62(4):1122–31. https://doi.org/10.1002/hep.27923.

    CAS  Article  Google Scholar 

  46. 46.

    Abrami L, Dallavilla T, Sandoz P, Demir M, Kunz B, Savoglidis G, et al. Identification and dynamics of the human ZDHHC16-ZDHHC6 palmitoylation cascade. Elife. 2017. https://doi.org/10.7554/eLife.27826.

    Article  PubMed  PubMed Central  Google Scholar 

  47. 47.

    Shi W, Chen X, Wang F, Gao M, Yang Y, Du Z, et al. ZDHHC16 modulates FGF/ERK dependent proliferation of neural stem/progenitor cells in the zebrafish telencephalon. Dev Neurobiol. 2016;76(9):1014–28. https://doi.org/10.1002/dneu.22372.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Cao N, Li J, Rao Y, Liu H, Wu J, Li B, et al. A potential role for protein palmitoylation and zDHHC16 in DNA damage response. BMC Mol Biol. 2016;17(1):12. https://doi.org/10.1186/s12867-016-0065-9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

The authors would like to give our sincere appreciation to the reviewers for their helpful comments on this article and research groups for the TCGA, which provided data for this collection. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Author information

Affiliations

Authors

Contributions

HW designed the overall study and revised the paper, XQH performed public data interpretation, XH drafted manuscript. DRS supervised the experiments. LNJ, MRQ and QZX participated in data collection, SYN, WZJ and WD contributed to data analysis, WJC and ZJX participated in the molecular biology experiments. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Wen Huang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

: Figure S1. Overall research design. Flow-process diagram presenting the process of comprehensive analysis. Figure S2. (A)The upset plot of gene interactions among the seven types of AS events in TCGA LIHC cohort. (B) The upset plot of gene interactions among the seven types of survival relevant AS events. Figure S3. LASSO coefficient of survival relevant AS events. (A) AA. (B) AD. (C) AP. (D) AT. (E) ES. (F) ME. (G) RI. Figure S4. A graph of the error rate of cross-validation. (A) AA. (B) AD. (C) AP. (D) AT. (E) ES. (F) ME. (G) RI. Figure S5. (A) Heatmap of the AA events PSI value in HCC. The color from red to blue shows a trend from high expression to low expression. (B) Distribution of AA prognostic signature risk score. (C) The survival status and duration of HCC patients in AA prognostic signature. (D) Heatmap of the AD events PSI value in HCC. The color from red to blue shows a trend from high expression to low expression. (E) Distribution of AD prognostic signature risk score. (F) The survival status and duration of HCC patients in AD prognostic signature. Figure S6. (A) Heatmap of the AP events PSI value in HCC. The color from red to blue shows a trend from high expression to low expression. (B) Distribution of AP prognostic signature risk score. (C) The survival status and duration of HCC patients in AP prognostic signature. (D) Heatmap of the AT events PSI value in HCC. The color from red to blue shows a trend from high expression to low expression. (E) Distribution of AT prognostic signature risk score. (F) The survival status and duration of HCC patients in AT prognostic signature. Figure S7. (A) Heatmap of the ES events PSI value in HCC. The color from red to blue shows a trend from high expression to low expression. (B) Distribution of ES prognostic signature risk score. (C) The survival status and duration of HCC patients in ES prognostic signature. (D) Heatmap of the ME events PSI value in HCC. The color from red to blue shows a trend from high expression to low expression. (E) Distribution of ME prognostic signature risk score. (F) The survival status and duration of HCC patients in ME prognostic signature. Figure S8. (A) Heatmap of the RI events PSI value in HCC. The color from red to blue shows a trend from high expression to low expression. (B) Distribution of RI prognostic signature risk score. (C) The survival status and duration of HCC patients in RI prognostic signature. Figure S9. (A) Kaplan–Meier curve presenting survival in AA prognostic signature. (B) ROC analysis of the risk scores in AA prognostic signature. (C) Kaplan–Meier curve presenting survival in AD prognostic signature. (D) ROC analysis of the risk scores in AD prognostic signature. (E) Kaplan–Meier curve presenting survival in AP prognostic signature. (F) ROC analysis of the risk scores in AP prognostic signature. (G) Kaplan–Meier curve presenting survival in AT prognostic signature. (H) ROC analysis of the risk scores in AT prognostic signature. Figure S10. (A) Kaplan–Meier curve presenting survival in ES prognostic signature. (B) ROC analysis of the risk scores in ES prognostic signature. (C) Kaplan–Meier curve presenting survival in ME prognostic signature. (D) ROC analysis of the risk scores in ME prognostic signature. (E) Kaplan–Meier curve presenting survival in RI prognostic signature. (F) ROC analysis of the risk scores in RI prognostic signature. Figure S11. (A) Univariate Cox regression analyses in AA prognostic signature. (B) Univariate Cox regression analyses in AD prognostic signature. (C) Univariate Cox regression analyses in AP prognostic signature. (D) Univariate Cox regression analyses in AT prognostic signature. (E) Multivariate Cox regression analyses in AA prognostic signature. (F) Multivariate Cox regression analyses in AD prognostic signature. (G) Multivariate Cox regression analyses in AP prognostic signature. (H) Multivariate Cox regression analyses in AT prognostic signature. (I) Univariate Cox regression analyses in ES prognostic signature. (J) Univariate Cox regression analyses in ME prognostic signature. (K) Univariate Cox regression analyses in RI prognostic signature. (L) Multivariate Cox regression analyses in ES prognostic signature. (M) Multivariate Cox regression analyses in ME prognostic signature. (N) Multivariate Cox regression analyses in RI prognostic signature. Figure S12. Kaplan–Meier survival analysis for multiple HCC subgroups according to the ALL signature stratified by clinical variables. (A-B) Stage. (C-D) T status. (E-F) Gender. (G-H) Tumor grade. (I-J) Age. (K) N status. (L) M status. Figure S13. (A) Comparison of the ESTIMATE score between risk score low/high groups. (B) Comparison of the immune score between risk score low/high groups. (C) Comparison of the tumor purity between risk score low/high groups.

Additional file 2

: Table S1. Univariate Cox regression analysis between AS events and survial. Table S2. Formula prognostic signature of HCC. Table S3. The list of splicing factor (SF) genes. Table S4. Results of correlation between AS events and SFs.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Xu, Q., Xu, H., Deng, R. et al. Immunological significance of prognostic alternative splicing signature in hepatocellular carcinoma. Cancer Cell Int 21, 190 (2021). https://doi.org/10.1186/s12935-021-01894-z

Download citation

Keywords

  • Hepatocellular carcinoma
  • Alternative splicing
  • Tumor immune microenvironment
  • Prognosis
  • Immunotherapy
\