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 (Figure S1A). These results showed that ES was the predominant splicing pattern meanwhile AD marked as the least frequent.
Identification of the survival-relevant AS events
With the help of Univariate Cox regression analysis, we identified a total of 3294 AS events which were significantly related with survival (P<0.05). The detailed description was recorded in Table S1. The intersecting sets of genes and splicing subtypes were delineated in the UpSet plot (Figure S1B). Among these subtypes of AS events, ES was the predominant pattern. The volcano map was generated to display the AS events (Figure 1A). The first 20 significant survival-related AS events from the seven subtypes were summarized in Figures 1B-H.
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 Figures S2A-S2G, S3A-S3G and 2A-2B. Then, we performed multivariate Cox analysis 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 (Figures S4A, S4D, S5A, S5D, S6A, S6D, S7A and 2C). The allocations of risk score (Figures S4B, S4E, S5B, S5E, S6B, S6E, S7B and 2D) and dot pot of survival status (Figures S4C, S4F, S5C, S5F, S6C, S6F, S7C 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 (Figures S8A, S8C, S8E, S8G, S9A, S9C, S9E and 2F; all P<0.05). To assess the prognostic value of risk signatures in HCC cohort, we performed ROC curve analysis. Area under curves of risk scores at 1-, 2- and 3‐year survival times were all more than 0.7, suggesting great sensitivity and specificity of their survival predictive ability (Figures S8B, S8D, S8F, S8H, S9B, S9D, S9F and 2G). Besides, results of univariate Cox model (Figures S10A, S10B, S10C, S10D, S10I, S10J, S10K and 2H) and multivariate Cox regression analysis (Figures S10E, S10F, S10G, S10H, S10L, S10M, S10N and 2I), suggesting risk scores could act an independent indicator in HCC.
We performed a stratification analysis 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 (Figures S11A, S11B). Similarly, ALL prognostic signature presented excellent prognostic prediction performance for patients in T1-2 or T3-4 status (Figures S11C and S11D), patients male or female gendered (Figures S11E, S11F), patients in 1-2 or 3-4 tumor grade (Figures S11G, S11H), patients aged <=65 years or >65 years (Figures S11I, S11J), patients in N0 status (Figures S11K), and patients in M0 status (Figures S11L). These results suggested that it can be an outstanding predictor independent from clinical parameters in patients with HCC.
Confirmation of Prognostic Value of ALL prognostic signature
We compared risk score among different subtypes according to clinical variables. The risk score increased significantly with advanced tumor grade (most P<0.05, Figures 3A), advanced clinicopathological stage (most P<0.05, Figures 3B) and advanced T category (most P<0.05, Figures 3C). To explore whether ALL prognostic signature was the best prognostic indicator among various clinical characteristics, age, gender, clinicopathological stage, tumor grade, T status, N status and M status were extracted as the candidate prognosis predictive factors. We integrated these clinical variables then conducted the AUC curve analysis for 1-, 2-, and 3-year OS and observed that risk signature obtained the most AUC (Figures 3D, 3E, and 3F). We generated a nomogram including risk score and clinicopathological stage to forecast prognosis of patients with HCC (Figure 3G). Age, gender and tumor grade were rejected out of the nomogram because of their AUCs were less than 0.6. Calibrate curves indicated powerful prognostic predictive ability of 1-, 2- and 3-year OS in our nomogram plot (Figure 3H-I).
Correlation of Risk Score with Tumor Immune Environment Characterization
To further examine whether risk score can act as immune indicators in both TCGA and ICGC HCC datasets, we performed 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). 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; Figures 4A-D). We found that low-risk patients obtained a higher stromal score compared with high-risk HCC patients (Figure 4F). However, there was no significant difference regarding to immune score, estimate score and tumor purity (Figures S12A, S12B and S12C). Subsequently, we distinguished distinction of the immune-related signatures between these two subgroups. Figure 4G and 4H showed that immune-related signature of each patient with corresponding immune scores in low-/high-risk groups. We observed that the infiltration of aDCs, Macrophages, Neutrophils, Tfh, Th1 cells, Th2 cells, Tregs, and some immune signatures like APC costimulation, T cells costimulation, check-point, HLA molecule expression level, inflammation-promoting, MCH class I expression and IFN response were significantly increased with decreased risk score (Figure 4I). The CIBERSORT algorithm results indicated that proportion of reseing Dendritic cells was negatively associated with risk score (Figure 4J). Above results indicated that ALL prognostic signature may provide a novel approach to elucidate the characteristics of immunity regulatory network in HCC.
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 [31-33]. Subsequently, we correlated six key immune checkpoint inhibitors genes (PDCD1, CD274, PDCD1LG2, CTLA‐4, HAVCR2, and IDO1) [27-29]. And we analyzed the correlation between ICB key targets and ALL prognostic signature to reveal the potential player of risk signature in the ICB treatment of HCC (Figure 5A). We found 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; Figures 5B-H). 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 (Figure 5G), suggesting ALL prognostic signature might act a nonnegligible role in the prediction of responsiveness to ICB treatment in patients with HCC
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. We compared ZDHHC16 expression level between normal tissues and tumor samples based on TCGA data. Relative to tumor tissues, ZDHHC16 expression level was lower in adjacent normal specimens (Figure 6A). Taking advantage of qRT-PCR, we determined ZDHHC16 expression level in two distinct HCC cell lines and human hepatic cell line. Consistent of results of online database, ZDHHC16 was upregulated in cancer cells relative to normal cell (Figure 6B). The expression level analysis among major pathological stages presented that ZDHHC16 expressed statistical significantly among different pathological stages (Figure 6C, F= 3.45 and P= 0.0168). And we found that the later tumor grade, the higher risk score (Figure 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 Figures 7E and 7F, lower ZDHHC16 expression level significantly suggested longer overall survival time (P = 0.0056) and longer disease-free survival time (P = 0.02). We compared 47 immune check blockade-associated genes expression levels between low-ZDHHC16 and high-ZDHHC16 groups and observed that 16 of 47 (i.e., PDCD1, CTLA4, etc.,) immune check blockade-associated genes expression levels were significantly dysregulated in between different subgroups (Figure 6G). Then we performed the correlation between the ZDHHC16 and ICB key targets adjusted by tumor purity by TIMER 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; Figures 6H-K), suggesting ZDHHC16 may exert a vital player in ICB treatment of HCC.
Development of the SF-AS regulatory Network
To elucidate the underlying mechanism of AS regulation, we generated a correlation network between the expression level of SFs and the PSI values of prognosis-related AS events. We identified 55 up-regulated AS events (yellow ellipses), 56 down-regulated AS events (blue ellipses) and 106 SFs (Figure 7). In our regulatory network, the top four most significant nodes were termed as the hub SFs or AS events (Table S), including one downregulated AS event (ACAA1-64022-ES), one upregulated AS events (SCP2-3045-ES) and two SFs (ISY1 and CLK2). As such, we speculated that these SFs served as pivotal regulators involving in the dysregulation of AS in HCC, further mediated tumor initiation and progression.