Our study included one hundred and twenty-five patients with hepatocellular carcinoma (HCC). Tumor puncture samples were collected from Qionglai Hospital of Traditional Chinese Medicine, and the whole exome sequencing was performed. The selected patients’ baseline demographic, clinical, and pathological data were collected and recorded during diagnosis, including sex, age, AFP, PIVKA-II, tumor size and number, tumor stage, treatment mode, survival status, and the latest follow-up information. The review date from the date of diagnosis to the date of death or the date of the last follow-up was collected for OS calculation, and the date from the first treatment to the time of the first recurrence was collected for PFS calculation. Exclusion criteria: when exploring whether patients with stage TNMI and TNMII relapse within one year, filter out patients with a follow-up period of less than one year; in the survival analysis of patients with stage TNMIII and TNMIV, filter out patients with an overall survival period of less than one month. Details were shown in Supplementary Table S1. The study was approved by the local ethics committee [No. 2021 (01)], and the patient’s informed consent was obtained following Chinese law.

Somatic mutation data were derived from TCGA (364 HCC samples) and ICGC (904 HCC samples, including LICA-CN, LICA-FR, and LIRI-JP) databases. Gene expression data (raw counts) were derived from TCGA (375 HCC samples, 50 normal tissue samples) and ICGC (LIRI-JP, 227 hepatocellular carcinoma samples). The original counts of gene expression data were normalized to CPM, and HCC samples with gene expression data, somatic mutation data, and complete clinical information were screened for subsequent analysis. Exclusion criteria: At the time of survival analysis, patients with overall survival of less than one month and a follow-up period of fewer than two months were excluded.

Whole exome sequencing was performed using the Agilent SureSelect Human All Exon v6 Kit (Agilent Technologies, Santa Clara, CA, United States), then sequenced on the Illumina HiSeq 2000 System (Illumina, San Diego, CA, United States) by Shanghai Personal Biotechnology Cp. Ltd. The 150bp paired reads were aligned to the human reference genome (GRCh38) by BWA (version 0.7.15). GATK (version 4.0) detected single nucleotide variants and indels. Copy number variation (CNV) was identified from the bam file by VarSeq (version 2.2.4, Golden Helix^®). All the variants were annotated by VarSeq with the public (gnomAD, 1,000 Genomes, DGV, ICGC, TCGA) and in-house databases. Variants passed quality control (Read depth ≥6, Genotype Qualities ≥20) were filtered. Variant filtering criteria: 1. Read depth filtering. A variant will be retained only if it meets the following conditions. (1) The total reads depth ≥10; (2) The alternative reads depth ≥5; (3) The alternative reads depth/total reads depth ≥0.1.2. Variant annotation result filtering. If a site meets the following conditions, it will be filtered. (1) The variants located upstream and downstream of the gene were filtered out; (2) The variants located in the intron region were filtered out; (3) Synonymous variants were filtered out; (4) The variants located on the transcript ablation region were filtered out; (5) The variant pathogenicity classified as “LOW” or “MODIFIER” were filtered out.

SPSS (Statistical Product Service Solutions, version 26.0) is used to analyze the correlation between molecular and clinical features and the statistics of clinical information. Rstudio (4.1.0) for data processing, analysis, and graphic visualization. GSEA is used for pathway enrichment analysis to determine the enrichment differences in biologically significant pathways among HCC patients in different groups. FDR <0.25, p < 0.05, and | NES | > 1 indicates that the results are meaningful. R package enrichplot for Gene Ontology (GO) and Jingdu Encyclopedia of Gene and Genome (KEGG) Analysis. R package Maftools is used to read and visualize MAF files, and R package GenVisR is used to draw a gene mutation waterfall map. R package edgeR and limma are used to identify differentially expressed genes (DEGs) among different groups of HCC patients. The cpm function in edgeR can convert the original data of gene expression into CPM format. R package survival and survminer are used for Kaplan-Meier survival analysis and graphical visualization. Evaluation of the predictive ability of R package timeROC and ggDCA in prognostic models. R package regplot for the Construction of Nomotu. R package ggplot2, pheatmap and ggpubr are mainly used for graphic visualization, including heat maps, histograms, volcanoes, etc. CIBERSORT algorithm is used to evaluate the relative proportion of 22 immune cells. The multi-algorithmic immune cell infiltration result file of the TCGA database is obtained at the website (http://timer.comp-genomics.org). The CMap database (http://clue.io/query) is used to screen small molecule drugs.

The correlation analysis between classified variables were processed by the chi-square or Fisher exact test. All statistical analysis was carried out using R software (version 4.1.0) or SPSS (version 26.0).

The median age of 125 patients with hepatocellular carcinoma was 53 years (range: 28–78 years), and the proportion of males was high (84.8% vs. 15.2%); Tumor staging using TNM classification was TNM I in 42.4% (n = 53), TNM II in 24.8% (n = 31), TNM III in 23.2% (n = 29), and TNM IV in 9.6% of the cases (n = 12); 75 resectable patients (60%) and 50 unresectable patients (40%) received curative surgery, immunotherapy, radiofrequency ablation (RFA), palliative treatment or combined treatment (Supplementary Table S1). The mean depth of target sequencing of the samples was 105X. The percentage of the target regions with a mean depth over 20X was 99.3%. Since we only get the WES data of tumor puncture samples from HCC patients, there may be errors in directly describing the gene mutations in our cohort. To better describe the mutations of the main driver genes of hepatocellular carcinoma (HCC), we screened the top 40 genes with mutation frequency from the somatic mutation data in the ICGC database and the TCGA database, including the four cohorts of ICGC-LICA-CN, ICGC-LICA-FR, ICGC-LIRI-JP and TCGC-LIHC, respectively, among which TTN, TP53, OBSCN, APOB, ADGRV1, XIRP2, PCLO, CSMD1, USH2A, LRP1B, FAT3, CSMD3, RYR2, and HMCN1 are repeated in four HCC queues (Figures 1A–D). In our cohort (N = 125), the frequencies of mutations in these genes were equally high, predominantly with missense mutations, and TTN was the gene with the highest mutation frequency; other genes with higher mutation frequencies were as follows: TP53 (30%), OBSCN (23%), APOB (22%), ADGRV1 (22%), etc. The mutation frequency was altered compared to the somatic mutation data for HCC in the TCGA and ICGC databases (Figure 1E). Analysis of tumor TNM stage showed that TP53 mutations and FAT3 mutations were distributed differently in different TNM stages (p < 0.05), and both were more common in patients with advanced HCC. Serum AFP level is a commonly used and vital index in diagnosing HCC and monitoring therapeutic efficacy. LRP1B mutations were mainly enriched in the normal serum AFP level (p = 0.044). We also found that APOB mutations were primarily distributed in HCC patients with tumors larger than 5 cm (p = 0.037) and XIRP2 mutations were primarily distributed in non-smoking patients (p = 0.002) (Figures 2A, C). By Fisher’s exact test, we found an evident coexistence between ADGRV1 mutations and LRP1B mutations, and between XIRP2 mutations and CSMD3 mutations (p < 0.05), there was significant mutual exclusion between the APOB mutation and the USH2A mutation (p < 0.05) (Figure 2B). In addition, we statistically analyzed mutations in main driver genes with other major clinical features and found no statistical significance (Supplementary Table S2).

To better distinguish our clinical patient samples, we focused on 1-year progression-free survival (PFS) of TNMI and TNMII HCC patients (Supplementary Table S3) and overall survival (OS) of TNMIII and TNMIV HCC patients (Supplementary Table S4). Through univariate COX regression analysis, we found that TP53 mutations, LRP1B mutations, sex, and serum abnormal prothrombin II (PIVKA-II) levels were closely related to recurrence within one year in TNMI and TNMII patients (p < 0.05) (Figure 2D). In contrast, TP53 mutations, LRP1B mutations, and serum AFP levels were closely related to death in TNMIII and IV patients (p < 0.05) (Figure 2G). In multivariate COX regression analysis, we found that TP53 mutations, LRP1B mutations, and serum abnormal prothrombin II (PIVKA-II) levels were independent factors of recurrence in HCC patients with TNMI and TNMII within one year (Figure 2E). TP53 mutations, LRP1B mutations, and serum alpha-fetoprotein (AFP) levels were independent influencing death factors in TNMIII and TNMIV HCC patients (Figure 2H). At the same time, the Kaplan-Meier survival analysis also validates our findings (Figures 2F, I). To further verify our conclusion, we performed a Kaplan-Meier survival analysis on the mutation data with complete clinical information in TCGA and ICGC database and also found that LRP1B mutations and TP53 mutations were significantly correlated with the prognosis of HCC patients (Figure 3A).

To further explore the biological characteristics related to TP53 mutation and LRP1B mutation, we downloaded RNA sequencing data from the TCGA database and screened out samples with RNA sequencing data and somatic mutation data. According to the mutation state of TP53 and LRP1B, the RNA sequencing data were divided into mutant group and wild group for follow-up analysis. Because many pathways are involved in tumorigenesis, the poor prognosis of HCC patients with TP53 or LRP1B mutations may be related to the abnormal regulation of many signal pathways in HCC. Through GSEA enrichment analysis, we found seven KEGG signal pathways, including DNA replication, base excision repair, homologous recombination, cell cycle, mismatch repair, splicing, and nucleotide excision repair, showed significant differential enrichment in TP53 or LRP1B mutant groups (Figures 3B, C). We then analyzed the infiltration of 22 kinds of immune cells in the tumor immune microenvironment of HCC patients. We found that memory B-cell, follicular helper T-cell, Treg cells, M0 macrophages, M1 macrophages, resident DC cells, and activated DC cells were significantly enriched in the TP53 mutation group, and CD4 resident memory T-cell were significantly improved in the TP53 non-mutation group (Figure 3D). Perhaps due to the small proportion of patients with LRP1B mutation in TCGA, we only found that activated NK cells were significantly enriched in the LRP1B mutation group, and the infiltration of other immune cells did not show a significant statistical difference (Figure 3E).

We created a detailed analysis flow chart (Figure 4A) to describe our analysis process more intuitively. In the TCGA database, 335 HCC patients have both RNA sequencing and somatic mutation data. Patients with HCC were divided into groups according to TP53 mutation status, including 96 TP53 mutant samples and 239 TP53 wild-type samples. Through differential gene expression analysis, we found that 809 genes were highly expressed in TP53 wild-type samples (logFC >1, p < 0.05), including 473 mRNA (Figure 4B). GO enrichment analysis showed that highly expressed genes were significantly enriched in the metabolic hormone process, muscle system process, collagen-containing, extracellular matrix apical part of cell, receptor ligand activity and signaling receptor activator activity (p < 0.05) (Figure 4C). KEGG enrichment analysis showed that highly expressed genes were significantly enriched in Neuroactive ligand-receptor interaction and Wnt signaling pathway pathways (Figure 4D). To ensure the integrity and rationality of clinical information, 167 samples of TP53 wild-type were included in the training set for constructing prognostic models. Through batch univariate COX regression analysis, it was found that 21 mRNA were closely related to the total survival time (OS) of patients with TP53 wild-type HCC (Figure 4E). Then, stepwise multivariate COX analysis, 11 mRNA were selected to construct the prognosis model of TP53 wild-type HCC patients to narrow the scope of gene screening further. The COX regression coefficient of related genes in each model is calculated based on COX multiple regression model, and the model risk score of each sample is defined as the product of the regression coefficient and the expression value of related genes in each model. According to the median risk score, the training set samples are divided into high-risk and low-risk groups for follow-up analysis. In addition, we further analyzed the distribution of the expression levels of 11 genes involved in constructing the model in different clinical characteristics, including hepatitis B virus infection, grade, stage, age, and sex (Figure 4F).

In the training set, the Kaplan-Meier survival analysis showed a significant difference in prognosis between the high-risk and low-risk groups. The median survival time of the low-risk group was significantly longer than that of the high-risk group (Figure 5A). Univariate and multivariate cox regression analysis was used to evaluate independent prognostic factors in patients with TP53 wild-type HCC. First, a univariate cox regression analysis was performed on all clinical variables and model risk scores. We found that age, sex, grade, and model risk score were closely related to the total survival time (OS) of patients with TP53 wild-type HCC (p < 0.05) (Figure 5B). Through further multivariate COX regression analysis, we found that the risk score of this model could be used as an independent prognostic factor for patients with TP53 wild-type HCC (p < 0.001) (Figure 5C). ROC analysis showed that the model could effectively predict the 1-year, 3-year, and 5-year survival of patients with TP53 wild-type HCC. The predictive ability was better than other clinical features (including age, sex, grade, stage, and hepatitis B virus). A higher AUC value indicates that the model has high sensitivity and specificity (Figure 5D). At the same time, DCA analysis also showed the superiority of the model’s predictive capacity for predicting 1-year, 3-year, and 5-year survival of patients with TP53 wild-type HCC (Figure 5E). The linkage analysis of risk factors showed that the low-risk and high-risk groups could be distinguished according to the median risk score of the model, and the death toll of the high-risk group was significantly higher than that of the low-risk group. And there are differences in the expression of genes involved in constructing the model in high and low-risk groups (Figure 5F). Next, we use TP53 wild-type HCC samples from the ICGC database (ICGC-JP) as a verification set further to verify the predictive ability of the prognostic model. Consistent with the above analysis method, we divided the validation set into high-risk and low-risk groups according to the model risk score. Kaplan-Meier survival analysis showed that the survival probability of the high-risk group was significantly lower than that of the low-risk group (Figure 5G). ROC analysis showed that the prognostic model could also effectively predict the 1-year, 3-year, and 5-year survival of patients with TP53 wild-type HCC (Figure 5H). The risk factor linkage analysis also draws the same conclusion as the training set (Figure 5I). As a supplement, we also evaluated the prognostic potential of the prognostic model in patients with mutant TP53. The risk score of each TP53 mutant patient in the TCGA cohort was calculated using the same algorithm and grouped according to the median risk score. Kaplan-Meier survival analysis showed that the prognostic model did not maintain a good predictive performance in patients with TP53 mutant HCC (Figure 6A). To sum up, the prognostic model is reliable in predicting the clinical prognosis of patients with TP53 wild-type HCC.

To explore which biological processes are involved in the poor prognosis of patients with TP53 wild-type HCC in the high-risk group, we performed a GSEA analysis. The results showed that the high-risk group was highly enriched in metabolic-related biological processes, such as oxidative phosphorylation and pentose phosphate pathway (Figure 6B). This suggests that the imbalance of metabolic processes may be a factor for poor prognosis in the high-risk group of patients with wild-type HCC of TP53. At present, immunotherapy plays an important role in tumor therapy, and the difference in immune-related characteristics of different HCC patients may affect the effect of immunotherapy. Below, we further analyze whether there are differences in immune-related characteristics between high and low-risk groups of TP53 wild-type HCC patients. The association analysis between joint immune checkpoints and high and low-risk groups showed that NRP1, CD160, and KIR3DL1 were highly expressed in the low-risk group, while CD276, TNFRSF14, LGALS9, TNFRSF4, and TNFRSF18 were highly expressed in the high-risk group (p < 0.05) (Figure 6C). However, the analysis of immune-related function only showed that Type_II_IFN_Reponse was enriched in the low-risk group (Figure 6D). We then used five algorithms: CIBERSORT, CIBERSORT−ABS, QUANTISEQ, MCPCOUNTER, and XCELL, to analyze the correlation between immune cell infiltration and high and low-risk groups. The results showed that T-cell regulatory (Tregs), Macrophage M0, Macrophage M2, Neutrophil, T-cell CD4⁺ (non-regulatory), Common lymphoid progenitor, T-cell NK, B-cell plasma, T-cell CD4⁺ Th1, T-cell CD4⁺ Th2 are enriched in the high-risk group, Monocyte, Mast cell activated, NK cell, Endothelial cell, Common myeloid progenitor, Granulocyte-monocyte progenitor, Hematopoietic stem cell are enriched in the low-risk group (Figure 6E). Considering the clinical utility of risk score in predicting OS in patients with TP53 wild-type HCC, we established a Nomograph including model risk score and main clinical features to predict 1-year, 3-year, and 5-year OS rates (Figure 6F). Next, we screened the small molecular compounds with potential therapeutic significance in high-risk patients with TP53 wild-type HCC. First, the differences in gene expression profiles between high and low-risk groups of TP53 wild-type HCC were analyzed, and the upregulated genes of the high-risk and low-risk groups were uploaded to the CMap database for comparison. We found that the expression profile of AM580 was similar to that of upregulated genes in the low-risk group (score = 99.01) but opposite to that in the high-risk group (score = −99.29) (Figure 6G). This suggests that AM580 has potential therapeutic significance for patients with TP53 wild-type high-risk HCC, and we provide a structural formula for AM580 (Figure 6H).