3.1 DEGs identification and prognostic model construction

The workflow for this study was shown in Fig. 1A. Firstly, based on the gene expression data of HCC from TCGA, DEGs involved in TCR signal pathway between HCC tumor and adjacent non-tumor tissues were identified. Overall, 92 DEGs were identified between tumor and adjacent non-tumor tissues, of which 52 genes were down-regulated in tumor compared with adjacent non-tumor tissues and 40 genes were up-regulated. These 92 identified DEGs were then subjected to univariate and multivariate Cox regression analysis. Eleven genes were identified to be potentially related to survival and might have clear effects on the prognosis of HCC OS (Fig. 1B). Using LASSO-Cox regression analysis, an 8-gene signature prognosis model was constructed (Fig. 1C). The risk score for prognosis prediction is calculated as follows:

Fig. 1

Model construction. A The workflow chart of this study. B Multi variate Cox regression analysis was used to identify DEGs. C LASSO regression was used to identify DEGs

Risk score = 0.370 × NRAS + (− 0.212) × PIK3R1 + (− 0.547) × AKT2 + (− 0.061) × CD247 + (− 0.086) × CD4 + 0.253 × CDC42 + 0.426 × CDK4 + (− 0.171) × PRKCQ.

Based on the median risk score, the HCC patients were assigned into low-risk and high-risk groups. The OS of the high-risk group was significantly lower than that in the low-risk group (P < 0.0001) (Fig. 2A).

Fig. 2

Model validation. A Using TCGA cohort as the training cohort. Validated at Nanfang Hospital (B) and ICGC cohort (C). D ROC analysis predicted OS using the risk score in Nanfang Hospital cohort. E RFS analysis

3.2 Well validation of the prognostic model across datasets

To validate the classification performance of the risk score signature, both the Nanfang Hospital and ICGC HCC cohorts were used as the external datasets. Likewise, patients were divided into two groups based on the median of the risk score. Obviously, the Kaplan–Meier survival curves also showed significant favorable OS in the groups with low-risk scores, and unfavorable OS in the groups with high-risk scores in both the ICGC and Nanfang Hospital HCC cohorts (P < 0.01) (Fig. 2B, C). The risk score value from Nanfang Hospital cohort was shown in Supplementary Fig. 1 and Supplementary Table 1. Meanwhile, according to the ROC analysis, the area under the curve (AUC) of the risk score for OS prediction in the Nanfang Hospital cohort was 0.728 at 1 year, 0.676 at 2 years, and 0.572 at 3 years (Fig. 2D). Additionally, we found out that there were more patients without recurrence in low-risk group than high-risk group (Fig. 2E). The recurrence free survival (RFS) in low-risk group was substantial longer than in high-risk group. Among the three different HCC cohorts, our prognostic model showed that it had the consistent OS predicting ability, particularly demonstrating more precise predictive competence for one year survival in the Nanfang Hospital HCC cohort.

3.3 The correlations between the prognostic model and HCC clinical features

To better elucidate the application of this model in clinical practice, we further analyzed the correlations between HCC clinical features and the prognostic model among the Nanfang Hospital cohort. We found out that tumor differentiation grade was lower in high-risk group than in low-risk group (Fig. 3A), which was consistent with the worse survival in high-risk group. On the contrary, the fibrinogen (Fbg) concentration, which is correlated with liver function and plays a vital role in coagulation cascade, was higher in high-risk group compared with low-risk group (Fig. 3B). Meanwhile, tumor ploidy (Fig. 3C) and tumor neoantigen burden (TNB) (Fig. 3D) were also higher in high-risk group than in low-risk group. Although there was no statistically significant difference of tumor mutation burden (TMB) and tumor purity between these two risk subgroups, the levels of TMB and tumor purity elevated slightly in high-risk group (Fig. 3E, F). Moreover, the multivariate cox regression analysis revealed that the low-risk scores, TNM (IV stage) and loss of heterogeneity in human leukocyte antigen (HLA) were independent prognostic factors for HCC patients (Supplemental Fig. 2A). The nomogram was then constructed (Supplementary Fig. 2B). In conclusion, by analyzing the prognostic model in connection with some common clinical indicators, the results were almost consistent with clinical characteristics, which further confirmed the applicability and effectiveness of our prognostic model for HCC patients.

Fig. 3

Clinical and tumor molecular features analysis in risk groups. A Tumor differentiation. B Fibrinogen concentration. C Tumor ploidy. D Tumor neoantigen burden (TNB). E Tumor mutation burden (TMB). F Tumor purity. G–K Top 5 gene mutation signatures in risk groups

3.4 Gene mutations in different risk groups

As for the demonstration of the molecular features among different risk groups, we compared the top 5 most common mutation genes (TP53, CTNNB1, MUC4, MUC16 and TNN) between high-risk and low-risk groups in Nanfang Hospital cohort (Supplementary Table 2). According to the Chi-square analysis, the TP53 mutation ratio in high-risk group was significantly higher than in low-risk group (69% vs. 8%), which was consistent with the prior study [16] (Fig. 3G). However, we found no statistical significance in CTNNB1, MUC4, MUC16 and TNN between these two risk groups (Fig. 3H–K).

3.5 The relationship between the prognostic model and TME

Given that the DEGs were analyzed based on TCGA data (Supplementary Fig. 3A), we further verified their expression in ICGC (Supplementary Fig. 3B) and Nanfang Hospital cohort (Fig. 4A). The expression levels of 8 genes in TGCA and IGCG cohort exhibited statistical significance between tumor and non-tumor tissues. Compared with non-tumor tissues, the expression levels of NRAS, CDC42 and CDK4 were higher in tumor tissues. In contrast, the expression levels of PIK3R1, CD247, CD4, AKT2 and PRKCQ were higher in non-tumor tissues compared with tumor. However, there was no difference in CDC42 expression between tumor and non-tumor tissues in Nanfang Hospital cohort. Furthermore, to clarify the characteristics of TME between different risk subgroups, we compared the expression levels of 8 genes respectively in tumor tissues among Nanfang Hospital cohort. Interestingly, the expression levels of NRAS, CDC42 and CDK4 were all elevated in high-risk group, while CD247 and PRKCQ were higher in low-risk group compared with that in high- risk group (Fig. 4B). In non-tumor tissues, only NRAS and AKT2 expressed higher level in high-risk group (Supplementary Fig. 4). Taken together, NRAS, CDC42 and CDK4 were considered to be associated with the poor prognosis of HCC patients.

Fig. 4

Expressions of 8 genes in prognostic models and immune features. A The expressions of 8 genes in prognostic model in Nanfang Hospital cohort. B The expression of 8 genes in tumor. C Immune Score in risk groups. D Immune cell infiltration in risk groups

Moreover, we analyzed the immune characteristics in these two subgroups. Firstly, we calculated the immune score in tumor and adjacent non-tumor tissues, respectively, then compared the T/N foldchange of it between the low-risk and high-risk groups. The T/N foldchange of immune score was significantly lower in the high-risk group than in the low-risk group (Fig. 4C). Furthermore, the features of 22 different types of immune cells in the different risk subgroups were also deliberated. The findings suggested that the T/N foldchange of naive B cells in the high-risk group was reduced compared with the low-risk group. Oppositely, the T/N foldchange of follicular helper T cells (Tfh) was significantly higher in the high-risk group (Fig. 4D). To conclude, those findings indicated that in high-risk group, immune infiltration in tumor was lower than that in adjacent non-tumor tissues.

3.6 TCR features in prognostic model

To further reveal the TCR features between the high-risk and low-risk group, we analyzed TCR clonotypes using MiXCR from RNA-seq data of Nanfang Hospital cohort. Firstly, the number of unique TCR-alfa chain (TRA) and TCR-beta chain (TRB) clonotypes in tumor tissues was significantly decreased in the high-risk group compared with the low-risk group (Fig. 5A), indicating that the TCR distribution was more abundant in tumor of low-risk group. Moreover, in high-risk group, the number of unique TRB clonotypes was higher in non-tumor tissues than tumor tissues. However, in non-tumor tissues, there was no difference in the number of unique TRA and TRB clonotypes between these two risk subgroups. Additionally, we also evaluated the Normalized Shannon diversity entropy (NSDE) value of TCR between high-risk and low-risk groups, but no differences were found (Supplementary Fig. 5A).

Fig. 5

TCR features in risk groups. A The number of unique TRA and TRB clonotypes. B The number of TRA expanded clonotypes. C The number of TRB expanded clonotypes. D–K The 10 top clonotypes in different subgroups

As we know, the expanding clones (ECs) may represent potential tumor-specific clonotypes. In this study, ECs were defined as those with a higher frequency than the average of unique clonotypes. Based on this, we analyzed the ECs of TRA and TRB in both tumor and non-tumor tissues. There was no difference in the average number of unique TRA or TRB clonotypes between high-risk and low-risk groups (Supplementary Fig. 5B, C). However, the number of TRA and TRB ECs in tumor was lower in high-risk group than in low-risk group, and especially in high-risk group, the number of TRB ECs was higher in non-tumor tissues than tumor tissues (Fig. 5B, C). All things considered, the TCR clonal expansion was more pronounced in low-risk group, and might oppositely be suppressed in tumor tissues of high-risk group.

A thorough analysis of the top 10 clonotypes was conducted among different subgroups. In the tumor group, we only found 1 TRA clonotype overlap (CAVMDSSYKLIF) between high-risk and low-risk group (Fig. 5D, E). However, in high-risk and low-risk subgroups of the non-tumor group, three identical and overlapping TRA clonotypes (CAMDSNYQLIW, CAVLDSNYQLIW and CAVRDSNYQLIW) were identified. These three clonotypes were likewise the top 3 clonotypes of non-tumor tissues in high-risk subgroup, and they were almost identical with the exception of the fourth amino acid (Fig. 5H, I). However, neither in tumor (Fig. 5F, G) nor non-tumor tissues (Fig. 5J, K) did we observe any overlap of the top 10 TRB clonotypes between high and low-risk groups. These top 10 TRA and TRB V/J gene ranks were shown in Supplementary Fig. 6 and 7, respectively.

3.7 BCR features in prognostic model

Similarly, the features of BCR immunoglobulin heavy chains (IgH) in different risk groups were also analyzed. In tumor, the number of unique BCR IgH clonotypes was significantly lower in high-risk group compared with low-risk group, which indicated there was less BCR IgH in tumor in high-risk group (Fig. 6A). At the same time, we observed that in both risk subgroups, the number of unique BCR IgH clonotypes in non-tumor tissues was higher than that in tumor (Fig. 6A). Likewise, even though there was no difference of NSDE value between high-risk and low-risk groups, the NSDE value of low-risk group was higher in non-tumor than tumor (Fig. 6B). Therefore, these findings demonstrated that BCR distribution trends varied among HCC risk subgroups. For further understanding the potential causes for the differences in regard to BCR, we analyzed BCR somatic hypermutation (SHM) and the expression of IgM/IgA related gene ‘JCHAIN’, which might be crucial to the function of B cells. We discovered that both their T/N foldchange were significantly or relatively lower in the high-risk group than that in the low-risk group (Fig. 6C, D). In other words, in the low-risk group, B cells in tumor experienced more antigen stimulation relative to the non-tumor than those in the high-risk group, which indicated the B cell function might be inhibited in the tumor of high-risk group.

Fig. 6

BCR features in risk groups. A The number of unique BCR clonotypes. B The NSDE value. C The BCR SHM foldchange comparison. D The JCHAIN expression. Gene-set enrichment analysis in high risk (E) and low risk subgroups (F)

Above all, these results indicated that in the high-risk group, the anti-tumor immune response for B cell function was suppressed in the TME of HCC.

3.8 Gene set enrichment analyses

To obtain the significant and targeted pathways involving related genes and clarify the potential molecular mechanism for the prognostic model, we conducted GSEA between different risk subgroups respectively in Nanfang Hospital cohort. It showed that there were different significant signaling pathway regulatory patterns in high-risk and low-risk groups. Firstly, the majority of the signaling pathways analyzed in this study were found to be upregulated in the high-risk group, with the exception of B cell receptor signaling pathway and Lysosome signaling pathway, which were shown to be downregulated (Fig. 6E). Nevertheless, in the low-risk group, only Protein digestion and abortion signaling pathway and GnRH secretion pathway were upregulated, whereas the others were under downregulation (Fig. 6F). Overall, in our prognostic model, the signaling pathway characteristics in different risk groups are almost completely different.