The transcriptome profiling RNA sequencing data [i.e., high-throughput sequencing data—fragments per kilobase of exon per million fragments mapped (HTSeq-FPKM)], the simple nucleotide variation data calculated by SomaticSniper software, and the clinical data were taken from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) data collection and were used to identify prognostic-related mutant genes. To validate the prognostic value of these mutant genes, simple nucleotide variation data and HCC clinical data were also downloaded from the International Cancer Genome Consortium Liver Cancer France (ICGC-LICA-FR) database. The patients in both databases with incomplete follow-up data were excluded. The ICGC and TCGA databases can be accessed without ethics approval as they are public databases; hence, we followed their guidelines for publication and data access policies in this study. A flow chart of the present study is shown in Figure 1 .

With the help of Perl software (version 3.6.1), gene mutation analyses of the ICGC and TCGA datasets were conducted using simple nucleotide variation data. Only mutations that caused a change in the sequence of amino acids were regarded as gene mutations and included in the subsequent analysis. Subsequently, according to the mutation statistical analysis, the waterfall plot of mutant genes in the above two datasets was mapped to visualize the genetic mutation in HCC samples, respectively, using the GenVisR software package.

To identify the prognostic-related mutant genes more accurately, the top 50 mutant genes of the ICGC and TCGA datasets were intersected. We conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses by using the clusterProfiler R package to explore the potential biological function of these intersected intersected mutant genes, including their cellular components (CCs), biological processes (BPs), signaling pathways, and molecular functions (MFs).

Tumor mutation burden (TMB) is defined as the number of mutations per 1 million bases. TMB analysis of the TCGA dataset was conducted using Perl software. In the present study, only a mutation that caused a change in the sequence of amino acids was regarded as a TMB and included in the subsequent analysis. The patients were categorized into the mutation group or the wild-type group based on the gene mutation analysis. The correlation analysis between TMB and intersected mutant genes (wild-type and mutation group) was conducted using the limma R package, and a box plot was mapped to visualize the difference using the ggpubr R package.

In the present study, the prognostic value of mutant genes was evaluated using the TCGA dataset and validated using the ICGC dataset. To identify the overall survival (OS)-related mutant genes, Kaplan–Meier (KM) analyses were conducted using the survival R package. Univariate and multivariate Cox analyses were conducted by utilizing the limma and survival R packages, and their clinicopathologic features and OS-related mutant genes were used to identify the HCC prognostic factors.

We conducted gene set enrichment analysis (GSEA) to explore the potential molecular mechanisms of OS-related mutant genes. To calculate immunocyte population fractions, the CIBERSORT method (13) was used, which can distinguish 22 phenotypes of immune cells sensitively and accurately. The samples with a p-value of < 0.05 were considered significant. According to CIBERSORT results, the correlation analysis among OS-related mutant genes and 22 phenotypes of immune cells was carried out using the limma and vioplot packages.

We developed a nomogram for predicting the 1-, 3-, and 5-year OS of patients with HCC by combining the prognostic-related factors selected from the univariate Cox regression analysis. We validated the predictive accuracy of the nomogram by comparing the observed actual probability with the calibration curve. The more consistent the reference line (black line) of the calibration curve (red line), the higher the accuracy of the nomogram.

A total of 24 HCC tissue samples and paired adjacent non-cancerous liver tissue samples were collected from patients who had undergone liver resection in The Second Affiliated Hospital of Guangzhou Medical University from May 2020 to December 2022. Before RNA and protein extraction, the samples were snap-frozen in liquid nitrogen and then stored at −80°C. This study was approved by the Ethics Committee of The Second Affiliated Hospital of Guangzhou Medical University (approval number 2023-ks-08).

The TRIzol™ reagent was used to isolate total RNA from the tissue samples. The SuperScript™ First-Strand Synthesis System was used to reverse-transcribe complement DNA (cDNA). Subsequently, real-time PCR (RT-PCR) was conducted, in accordance with the manufacturer’s instructions, by using an RT-PCR system (LightCycler^® 480 Instrument II; Roche Holding, Basel, Switzerland). The primers for the glyceraldehyde 3-phosphate dehydrogenase (GAPDH) gene were used as an internal loading control. The primers used were: DNAH5 (divergent primers), forward: 5′-GGACCTGGAGTTGCTGCTTGAC-3'; DNAH5, reverse: 5′-AGGCGTGCTGCTCATTTCTTCTAG-3'; GAPDH, forward: 5'-CTCCTAAACCATAGCACCATCG-3'; and GAPDH, reverse: 5′-CGTCAGTGCGGTCACAATC-3'. The 2^–ΔΔCT method was used to standardize all the values.

Radioimmunoprecipitation assay (RIPA) lysis buffer was used to extract proteins from the tissue samples, and a bicinchoninic acid (BCA) kit was used to quantify the proteins. Sodium dodecyl-sulfate polyacrylamide gel electrophoresis (SDS-PAGE) was carried out using a 7.5% gel to separate the proteins. The separated proteins were transferred onto a polyvinylidene difluoride (PVDF) membrane. Subsequently, the PVDF membrane was blocked in 5% skim milk for 2 hours at room temperature and incubated at 4°C overnight with a polyclonal rabbit anti-DNAH5 antibody (1:1,000; Invitrogen Life Technologies, Grand Island, NY, USA). The membranes were then incubated with horseradish peroxidase (HRP)-conjugated goat anti-rabbit secondary antibody (1:5,000; Abcam, Cambridge, UK) for 1.5 hours and visualized using the Immobilon™ Western Chemiluminescent HRP Substrate (Millipore, USA).

All bioinformatic and statistical analyses were conducted using R software (version 4.2.2) and Perl software (version 3.6.1). The genes with the top 50 mutation frequencies in the TGCA and IGGC databases were extracted using Perl. The R package “GenVisR” was used to visualize the mutations of these genes (14). These genes were intersected to obtain the genes with a high mutation frequency in both databases using the R package “venn.” The correlation analysis between the TMB and intersected mutant genes was assessed and visualized via the limma and ggpubr R packages, respectively. The GSEA analysis was conducted using the DNAH5 mutation and expression matrix data in the GSEA software (v. 4.1.0) (15). The normalized enrichment score (NES) was calculated by setting the permutation values to 1,000, and a false-discovery rate (FDR) p-value < 0.05 was adopted to determine the evident enrichment pathways. CIBERSORT is a computational method that can assess the proportion of 22 immunocyte immune cells in tumor tissue based on transcriptome data (13). The matrix data of the immune cell proportion for each tumor sample were obtained using the CIBERSORT deconvolution algorithm, setting the filter condition to a p-value < 0.05. The matrix data visualization was conducted using the R package “corrplot”. The TCGA samples were assigned to either the wild-type group or the mutation group based on their DNAH5 status. The difference analysis of immunocyte infiltration within the two groups was carried out using the R package “limma” and visualized using the R package “vioplot.” The survival curves were analyzed via KM survival analysis and evaluated using the log-rank test. The identification of prognosis risk factors was conducted through a survival analysis of the clinical characteristics of patients, including their age, gender, grade, stage, TMB, BCLC stage, and DNAH5 classification via univariate and multivariate Cox regression analyses. The correlation between the mutant genes and TMB was evaluated using the Mann–Whitney U-test. For all comparisons, a two-tailed p-value < 0.05 was considered statistically significant.

The gene mutation analysis showed that 10,981 and 11,098 mutant genes were identified in the ICGC and TCGA datasets. CTNNB1, TNN, and TP53 were the top three mutant genes in both the ICGC and TCGA datasets. The waterfall plots of the top 50 mutant genes in the ICGC and TCGA datasets are shown in Figures 2A, B , respectively. A total of 28 mutant genes were identified in both the ICGC and TCGA datasets by way of intersected mutant gene analysis, which is shown in Figure 2C . The GO analysis showed that “maintenance of location,” “myofibril,” and “actin filament binding” were the most common functions of these mutant genes in terms of their BPs, CCs, and MFs, respectively ( Figure 3 ). According to the KEGG analysis, these mutant genes were involved in the pathways of multiple cancers, including HCC, gastric cancer, breast cancer, and thyroid cancer ( Figure 3 ).

The correlation analysis between TMB and intersected mutant genes showed that 23 mutant genes had a higher TMB, and only five mutant genes had no correlation with TMB, namely, albumin (ALB), apolipoprotein B (APOB), CUB and sushi domain-containing protein 3 (CSMD3), spectrin repeat-containing nuclear envelope protein 2 (SYNE2), and spectrin alpha, erythrocytic 1 (SPTA1); this is shown in Figure 4 . To identify OS-related mutant genes, KM analysis of the above 23 TMB-related mutant genes was carried out. According to the results of the KM analysis, only patients with the DNAH5 mutation had poorer OS than those without the DNAH5 mutation in both the ICGC and TCGA datasets ( Figure 5A ). Univariate and multivariate Cox analyses of clinicopathologic features (including age, gender, tumor grade, pathological stage, and Barcelona Clinic Liver Cancer stage) and the mutant gene DNAH5 were conducted to identify independent prognostic factors. Finally, DNAH5 and pathological stage were independent prognostic factors in the TCGA dataset ( Figure 5B ), whereas in the ICGC dataset, DNAH5 and Barcelona Clinic Liver Cancer (BCLC) stage were independent prognostic factors ( Figure 5C ).

The GSEA (c2.cp.kegg.v7.symbol.gmt) revealed that the main biological functions of mutant gene DNAH5 were related to cancer and that all 15 pathways were highly expressed in the wild-type group compared with the mutant group, namely, apoptosis, the B-cell receptor signaling pathway, complement and coagulation cascades, cell cycle, DNA replication, colorectal cancer, mitogen-activated protein kinase (MAPK) signaling pathway, vascular endothelial growth factor (VEGF) signaling pathway, mammalian target of rapamycin (mTOR) signaling pathway, wingless-related integration site (Wnt) signaling pathway, Notch signaling pathway, transforming growth factor beta (TGF-β) signaling pathway, p53 signaling pathway, pathways in cancer, and the peroxisome proliferator-activated receptor (PPAR) signaling pathway ( Figure 6A ). The highly expressed pathways of the mutant group contributed to promoting cancer development and progression, whereas the low-expressed pathways of the mutant group contributed to inhibiting cancer development and progression. The CIBERSORT method was conducted to explore the relationship between the mutant gene DNAH5 and immune cells, as shown in Figures 6B, C . As a result, the mutant gene DNAH5 was negatively related to naive CD4 T cells, monocytes, activated mast cells, and activated dendritic cells compared with those without a mutation of the DNAH5 gene ( Figure 6D ).

We conducted paired differential expression analysis to explore the differential expression of DNAH5 between HCC samples and non-tumor liver samples in the TCGA and ICGC datasets. A total of 50 paired tissue samples were included in the TCGA dataset and 198 paired tissue samples were included in the ICGC dataset. The results show that in the ICGC dataset, the level of mRNA expression of DNAH5 in HCC samples was significantly higher than in the non-tumor samples, but no statistical difference was found in the samples from the TCGA dataset ( Figures 7A, B ).

Hence, to further explore the expression profile of DNAH5 in HCC, the mRNA and protein expression profiles of DNAH5 in HCC were validated in 24 pairs and 6 pairs of HCC tissue samples and adjacent non-cancerous liver tissue samples, respectively. As a result, the levels of mRNA and protein expression of DNAH5 were significantly higher in the tissue samples of patients with HCC ( Figures 7C–E ). The clinical features of those patients are shown in Supplementary Table 1 .

We integrated the mutant gene DNAH5 with other characteristics of HCC, including pathological stage and TMB, to develop a nomogram capable of predicting the 1-, 3-, and 5-year OS of HCC patients ( Figure 8A ). The calibration plots of this nomogram indicated that our nomogram performed well ( Figures 8B–D ).