Identification and Validation in a Novel Classification of Helicase Patterns for the Prediction of Tumor Proliferation and Prognosis

Introduction

The most common and lethal type of primary liver cancer is hepatocellular carcinoma (HCC), which accounts for approximately 75% to 90% of all cases and is followed by cholangiocarcinoma.^1–3 HCC is now the third leading cause of cancer-caused death worldwide, placing a heavy economic burden on healthcare systems and demonstrating that HCC is a major health problem.⁴ Although classical treatment strategies based on clinical risk stratification are currently widely used in the clinic, the effect of these strategies on improving the prognosis of HCC patients was recently challenged.^5,6 With the development of precision medicine based on molecular diagnosis and treatment, more and more attention is being focused on providing individualized treatment strategies for patients, especially on the discovery of new stratification methods and biomarkers.^7–9

Recent studies have shown that there are close direct interactions between helicases and several nuclear proteins that maintain chromosome stability and are potentially important in cancer susceptibility.¹⁰ As a key enzyme that determines cell division as well as nuclear transcription, helicases are easily associated with a unique and crucial role in DNA damage response, DNA repair and maintenance of genomic stability.¹¹ Considering that selective inactivation of DNA repair pathways can affect the efficacy of antitumor therapy, selective killing of cancer cells through inactivation of the resolvase-dependent biological pathway using means that regulate helicase expression or function is a feasible strategy.^12,13 Overexpression of DDX1 has been proved to be an independent prognostic indicator of early recurrence of breast cancer.¹⁴ In addition, it has been shown that high levels of DDX3 are associated with poor prognosis in gallbladder and breast cancers.^15,16 In contrast, there is little knowledge of the biological processes involved in DDX56 in cancer species. Only recently, the expression of DDX56 has been proved to be related to lymphatic invasion and distant metastasis in colorectal cancer,¹⁷ and it is frequently up-regulated in tumor tissue and osteosarcoma cell lines.¹⁸ DDX56 may be related to the early recurrence of squamous cell lung cancer, involving transduction of the Wnt signal pathway,¹⁹ which indicates the potential carcinogenic effect of DDX56 in cancer.

In this study, we aim to develop a helicase-related stratification strategy to predict the immune microenvironment, prognosis and response to immunotherapy in HCC and to further search for valuable diagnostic and prognostic markers. In the future, this work could help physicians make important judgments about diagnosis and treatment.

Materials and Methods

Reagents

All reagents and antibodies involved in this study are listed in the Supplementary Materials (Supplementary Table S1). The concentrations of antibodies used are based on those used in previous studies or as recommended by the relevant manufacturer. Specific experimental procedures were included in the Supplementary Materials.

The Research Flowchart

The research flowchart demonstrated our study design (Figure 1).

Figure 1 The flow chart showing the scheme of our study for hepatocellular carcinoma.

Clinical Sample Collection and Immunohistochemical (IHC) Staining

This study was conducted in accordance with the Declaration of Helsinki and approved by the ethics committee of the Zhangjiagang TCM Hospital Affiliated to Nanjing University of Chinese Medicine (Ethics Approval No.2022-03-14). Written informed consent was obtained from all patients. Inclusion/exclusion criteria for participants were provided in the Supplementary Material. Immunohistochemical Staining procedures and H-scoring methods have been previously described.^20,21 The detailed description was also included in Supplementary Materials.

Cell Culture and Western Blotting

Five human HCC cell lines were used in this study. Those human hepatocarcinoma cell lines (SNU-182/SNU-387/HepG2/Huh7/SK-Hep-1) came from a single supplier, Wuhan Pu-nuo-sai Life Technology Co. Ltd. (Wuhan, China). They were cultivated in Dulbecco’s Modified Eagle’s Medium (DMEM) supplemented with 10% fetal bovine serum (FBS). All the cells were maintained in a humidified 5% CO₂ atmosphere at 37°C. The Western blotting protocol was consistent with previous studies.²² The detailed description is shown in Supplementary Materials.

Public Datasets

The RNA-seq transcriptome information and matching clinical data of 374 Liver hepatocellular carcinoma (LIHC) patients and 50 normal control samples were acquired from TCGA portal website.²³ To ensure the reproducibility of the analysis results, independent datasets were downloaded from the Gene Expression Omnibus (GEO) database,²⁴ where GSE14520 (n = 220) was used for independent clustering, and GSE36376 (n = 433) and GSE102079 (n = 257) datasets were used for validation of the differential expression analysis results. The single-cell RNA data included in this study were obtained based on a previous HCC study (GSE125449).²⁵

Consensus Clustering

The “ConcensusClusterPlus” R package was conducted to perform consensus clustering to stratify HCC patients based on helicases.²⁶ The ideal cluster numbers between k = 2–6 was then assessed, and this procedure was repeated 1000 times to guarantee stability.

Identification of Differentially Expressed Genes in HCC

The differentially expressed genes (DEGs) were obtained by “Limma” package.²⁷ DEGs with an absolute log2 fold change (FC) >2 and an adjusted P value of <0.05 were considered with statistical significance, and then were used for subsequent analysis.

Enrichment Analysis

The “clusterProfiler” package was used for the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis.^28,29 In addition, the gene set enrichment analysis (GSEA)-based differential signaling pathways were also obtained by “Limma” package,³⁰ and the Hallmark gene set data resource derived from the Molecular Signatures Database-MsigDB (www.broadinstitute.org/gsea/msigdb) was downloaded as a background set. P < 0.05 was defined as significant enrichment.

Immune Analysis

CIBERSORT was used to quantify the level of infiltration of different immune cells in HCC, and tumor immune dysfunction and exclusion (TIDE) score was calculated to predict the potential immune checkpoint (ICB) blockade response in patients with HCC.³¹

Construction of Co-Expression Networks

The “multiscale embedded gene co-expression network analysis (MEGENA)” package in R software was conducted to identify a co-expression network.³² As an recently developed co-expression module analysis method, MEGENA offers unique advantages over traditional co-expression analysis tools in efficiently constructing large-scale co-expression plane filtering networks and preserving gene interactions.³² Specifically, fast planar filtered network (PFN) construction is the first step of MEGENA analysis, then significant gene pairs in PFN were obtained computationally, and finally the constructed PFNs are aggregated into multiscale clustering analysis (MCA) for subsequent analysis.³³

We extracted the largest gene module from the co-expression network and converted them into a format readable by Cytoscape software for final analysis as well as visualization.³⁴ Degree values were calculated and used to rank the genes in the module to identify potential hub genes.

Single Cell Analysis

Seurat (Version 3) (http://satijalab.org/seurat/) was used to process raw gene expression matrices. The filtering thresholds are set as follows: 1) Genes expressed in at least three cells. 2) Cells contain 500–6000 features. 3) Mitochondrial genes in cells are less than 20%. The highly variable genes (HVGs) were identified with “FindVariableFeatures” function. Principal component analysis (PCA) was conducted to determine the principal components suitable for subsequent analysis,³⁵ and uniform manifold approximation and projection (UMAP) was carried out to visualize the clusters.³⁶ The “FindAllMarkers” function was applied to find marker genes for each cell cluster. Clusters were annotated to known cell types by “singleR” package according to well-known marker genes.

GSCA Database

GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/) was used to analyze expression, DNA methylation, CNV (copy number variation), and drug sensitivity about helicases. The degree of helicase activation or inhibition of the cancer-related pathways was also analyzed.

Rusults

The Genetic Characteristics and Transcriptional Variations of 112 Helicases

Basic information and list of all helicases was obtained from the Human Genome Organization (HUGO) portal,³⁷ and 112 helicases mainly include DNA helicase and RNA helicase. In order to observe the genetic variation of helicase molecules in cancer, we showed somatic mutations of 10 molecules with the highest mutation frequency in cancer tissues according to TCGA pan-cancer data. Among 4798 samples, 2125 samples carry helicases molecular mutation, with a frequency of 44.29%. In all helicases, ATRX showed the highest mutation frequency (Figure 2A). Genetic variation was further confirmed to be an important factor that affects the expression of the helicases. For most helicases, there is a positive correlation between CNV and mRNA expression levels, while methylation and mRNA expression levels showed opposite trends (Figure 2B and C). Besides, we found that different helicases have very different prognostic significance in different cancer types (Figure S1). In the TCGA-LIHC dataset, we observed extensive and significant co-mutations between helicase molecules (Figure 2D). Among 363 samples, 176 samples had mutated helicase molecules, at a frequency of 48.48%. Among the 176 samples, DHX9 had the highest mutation frequency of up to 3% (Figure 2E), mainly missense mutation.

Figure 2 Expression variation of helicase molecules. (A) The waterfall diagram shows the somatic mutations of the 10 helicases with the highest mutation frequency in pan-cancer. Here, 44.29% is the proportion of 2125 samples with at least one mutation of the top 10 genes among 4798 samples with at least one mutation of 112 helicases. The percentage figure of each line on the right of the picture is the number of samples with the corresponding gene mutation divided by 4798 samples with at least one mutation among the 112 helicases. (B) The bubble chart shows the correlation between CNV of the 112 helicases and mRNA expression level. Red indicates positive correlation; blue indicates negative correlation. The deeper color indicates a larger correlation index. The bubble size indicates the FDR. (C) The bubble chart shows the correlation between methylation of the 112 helicases and mRNA expression. Red shows a positive correlation and blue shows a negative correlation. The darker color indicates a larger correlation index. Bubble size indicates the FDR. (D) Mutation characteristics of the 112 helicases in 374 patients with bladder cancer in the TCGA-LIHC cohort; green indicates co-mutation, the asterisk indicates P-value (P < 0.05, *P < 0.1). (E) Mutation frequency of 112 helicases in 374 patients in the TCGA-LIHC cohort. Each column represents an individual patient. The small figure above shows the TMB, the number on the right shows the mutation frequency of each regulator, and the figure on the right shows the proportion of each variant.

Two Different Helicases Level Patterns Identified by Unsupervised Learning

To fully understand the role of helicase in LIHC, we used consensus clustering to identify two unique modification patterns of LIHC associated with helicase. To fully understand the role of helicase in LIHC, we used consensus clustering to identify two unique modification patterns of LIHC associated with helicases, called Cluster 1 (C1) and Cluster (C2) (Figure 3A–C). Principal component analysis confirmed that these two clusters could be accurately distinguished by helicase expression profiles (Figure 3D). Overall, C2 showed high helicase expression levels, indicating a high-helicase subtype. In contrast, C1 showed low expression levels, indicating a low-helicase subtype (Figure 3E). In addition, survival analysis showed that the 2 helicase-based subtypes presented different clinical outcomes, with C1 having a poor prognosis and C2 having a favorable clinical outcome (Figure 3F, OS, P = 0.0205; Figure 3G, DSS, P = 0.0162; Figure 3H, disease-free survival (DFS), P = 0.134). The analysis procedure for GSE14520 was identical to TCGA-LIHC cohort (Figure S2). Interestingly, we observed a higher proportion of patients with advanced stage in C1 than in C2, especially in terms of T stage and pathological stage (P < 0.05, Figure 3I). We subsequently compare the difference in mutation frequency between the two subtypes and noted that TP53 was mutated much more frequently in C1 than in C2 (Figure 3J), implying that C1 has higher genomic instability. In addition, we found more mutational co-occurrence in C1, such as RYR1 and RYR2 mutations (Figure 3K), while in C2, most mutational events were relatively independent and did not interact with each other (Figure 3L).

Figure 3 (A) Heatmap depicts consensus clustering solution (k = 2) for 112 helicases in 374 LIHC samples; (B) Delta area curve of consensus clustering indicates the relative change in area under the cumulative distribution function (CDF) curve for k = 2–6. (C) The consensus score matrix of all samples when k = 2–6. (D) The PCA distribution of TCGA-LIHC cases by expression profile of helicases. Each point represents a single sample; different colors and shapes represent the C1 and C2 subtypes, respectively. (E) Heatmap of 112 helicases expressions in different subtypes. Red represents high expression and blue represents low expression. (F–H) Kaplan–Meier curves of OS (F), DSS (G), DFS (H) in C1 and C2 (TCGA-LIHC, n = 374). Log rank test was performed. (I) Relationship between two subtypes and the clinicopathological parameters, including gender, T, N, and M stage, pathological stage, and tumor grade (TCGA-LIHC, n = 374, *P < 0.05). (J) Oncoplot displaying mutation profile of a union set of the top 20 genes with highest mutation frequency in C1 and C2, respectively. (K and L) Graphs showing mutually exclusive and co-occurrence of mutational events in C1(K) and C2 (L), respectively (·P < 0.05, *P < 0.1). (M) Differential TIDE score between C1 and C2 (TCGA-LIHC, n = 374, **P < 0.1). Wilcoxon test was performed.

Finally, we analyzed the effect of helicases on tumor immune microenvironments. We use CIBERSORT to calculate the level of immune cell infiltration in LIHC samples using RNA-sequencing data. There was a higher quantity of M0 macrophage in C1 than in C2 (Figure S3A, P < 0.001), and the high infiltration of M0 macrophages was a factor in the poor prognosis of patients with LIHC (Figure S3B, OS, P = 0.02). Notably, all common ICPs (immune checkpoints) were significantly higher in C1 than in C2 (Figure S3C, P < 0.05), while the TIDE algorithm showed that patients in C1 may have worse immunotherapy outcomes (Figure 3M, P = 0.0051).

Identification of Differentially Expressed Genes (DEGs) and Signal Pathways in Different Helicases Subtypes

As the two subtypes exhibited different clinical outcomes and pathological features, we here identify key DEGs and signaling pathways in each subtype to understand the underlying regulatory mechanisms of the helicase phenotype. We found an aggregate of 921 DEGs between two clusters (Figures 4A), and the upregulated genes in C1 were mainly enriched in activities associated with fatty acid metabolism, spermatogenesis, hedgehog signaling, pancreas beta cells, KRAS signaling, inflammatory response (Figure 4B). To identify key gene modules, we constructed a MEGENA network using these DEGs (Figure 4C). The largest module C1_4 consisting of 358 genes was shown in Figure 4D. Enrichment analysis shows that this module is closely associated with nuclear division as well as genomic stability (Figure 4E). We subsequently conducted Matthews correlation coefficient (MCC) algorithm to identify 10 hub genes which had a consistent activation of the cell cycle (Figure 4F), and all hub genes significantly highly expressed in hepatocellular carcinoma tissues (Figure 4G, P < 0.001) and correlated with unfavorable prognosis (Figure S4, P < 0.05). In addition, we found that overexpression of most hub genes in the GDSC database was associated with higher drug sensitivity, except for KIF23 (Figure S5A). In the Cancer Therapeutics Response Portal (CTRP) database, overexpression of KIF18B, GINS1, CDCA5, KIFC1, KIF2C, KIF4A, and TOP2A implied insensitivity to Trametinib, a targeted drug that inhibits cell proliferation through MEK signaling (Figure S5B). Finally, we observed higher expression of 10 hub genes in C1 than in C2 (Figure 4H, P < 0.001), further suggesting that they are potential regulators of helicases subtypes differences.

Figure 4 Function analysis of helicase phenomenon. (A) Volcano map of differentially expressed genes between C1 and C2. Red indicates up-regulated genes in C1, blue indicates down-regulation. Data on the abscissa are differences in gene expression (log2 fold change); data on the ordinate represent the significance of these differences (−log10 padj) (TCGA-LIHC, n = 374). (B) Gene Set Enrichment Analysis (GSEA) of C1 and C2. (C) The DEGs between C1 and C2 in the co-expression network. Each node represents a gene module, with the larger nodes indicating a higher number of genes. (D) The MEGENA network showing the largest gene module. The degree value of the node was calculated by Cytoscape. The higher the degree value, the darker the color the larger the size. (E) GO and KEGG enrichment analysis of the largest gene modules. Each colored bar indicates a different biological process. (F) The correlation between the 10 hub genes in HCC and important cancer signaling pathways. The solid line represents activation and the dashed line represents inhibition. (G) Expression level of 10 hub genes between HCC tissues in the TCGA-LIHC cohort (n = 50) and paired control tissues (TCGA-LIHC, n = 50, ***P < 0.001). (H) Expression level of 10 hub genes between C1 and C2 (TCGA-LIHC, n = 374, ***P < 0.001).

Identification of Hub Genes Upregulation in Highly Proliferative Cells

Here, we sought to localize the hub gene at the single-cell level to investigate the biological mechanism of the helicase phenotype. The process of the single-cell analysis was demonstrated in Figures 5A and B with a Uniform Manifold Approximation and Projection (UMAP) and a t-Distributed Stochastic Neighbor Embedding (tSNE) that displayed the cell markers of each cell cluster. We then used singleR for automatic annotation of cell types (Figure 5C and D). The expression levels of CDCA5, KIF2C, KIFC1, KIF4A, HJURP, TOP2A, BUB1B, KIF18B, KIF23, GINS1, and MKI67 in a single-cell study are illustrated in Figure 5E and F. We observed that most hub genes are expressed in concert with Ki67, mainly on cells in the G2M phase (Figure 5G). Finally, in the TCGA-LIHC dataset, we further confirmed the significant positive correlation between MKI67 level and all hub gene levels (Figure 5H). These results suggested that the helicase phenotype was closely related to cell proliferation.

Figure 5 Single cell analysis of hub genes. (A) Heatmap depicting expressions of top 10 marker genes among 14 detected HCC cell clusters. (B) Cells were clustered into 14 types via UMAP and tSNE dimensionality reduction algorithm, each color represents a unique cluster. (C) SingleR was used to automatically annotate each cell, and the bar at the top of the graph represents the cluster and cell label, respectively. (D) The dot plot shows the situation after manual modification, with each color representing the annotated phenotype of each cluster. (E and F) tSNE distribution of all hub genes, including dot plot (E) and violin plot (F). (G) The cell cycle is annotated by marker genes, red for G1 phase, green for G2M phase and blue for S phase. (H) Heat map showing the significant correlation of all hub genes with MKI67. The top half of the graph represents the expression trend of MKI67 and the bottom half represents the expression trend of all hub genes, with blue representing low expression and red representing high expression (TCGA-LIHC, ***P < 0.001).

Identification of DDX56 Overexpression as a Poor Prognostic Factor for LIHC

In this part, LASSO regression of overall survival (OS), disease-free survival (DSS) and progression-free interval (PFI) was conducted based on 112 helicases, respectively, to identify robust markers for follow-up studies. Cross-validation was performed in ten rounds to prevent overfitting (Figure 6A–C). Sequential Cox regression (both univariate and multivariate) indicated that DHX58, DDX43, DDX56 were significantly correlated with LIHC prognosis (Figure 6D–I). Then, we constructed Kaplan–Meier survival curves to validate the prognostic value of DHX58, DDX43, DDX56 and the results showed that patients with higher DDX56 expression level had shorter PFI (Figure 6J–M, P < 0.05). In addition, we observed a trend of overexpression of DDX56 in almost all common tumor types (Figure 6N). These results suggested that DDX56 is a potential oncogene.

Figure 6 Identification of DDX56 as an independent prognosis factor. (A–C) least absolute shrinkage and selection operator model (LASSO) fitting based on OS (A), DFS (B), PFI (C). Each curve represents a gene. The profiles of coefficients were plotted versus log(λ). Vertical lines indicate the positions of all genes with coefficients greater than 0 determined by 10-fold cross-validation. λ was determined from 10-fold cross-validation. The x-axis represents log(λ); the y-axis represents binomial deviance. Optimal values calculated from minimum criteria and one standard error of the criteria are indicated by the dotted vertical lines. (D–I) Univariate (D–F) forest plot and multivariate (G–I) forest plot showing association between 6 candidate genes expression and OS (D and G), 8 candidate genes and DSS (E and H), 5 candidate genes and PFI (F and I) in HCC. (J–M) OS about DHX58, DSS about DDX43, PFI about DDX56, and DSS about DHX58 based on TCGA-LIHC (TCGA-LIHC, n = 374, Log rank test). (N) The mRNA expression of DDX56 between tumor and normal control tissues in pan-cancer was assessed from TCGA database. (TCGA-LIHC, n = 374; NS: not significant, *P < 0.05, ***P < 0.001).

Expression Analysis of DDX56 in LIHC

We first found that DDX56 was highly expressed in C1 with relatively poor prognosis (P = 1.7e-14, Figure 7A) and confirmed the overexpression of DDX56 in tumor tissue in paired TCGA samples as well as in other two independent datasets (p < 0.001, Figure 7B and C). We next examined the expression of DDX56 in several common hepatocellular carcinoma cell lines, and the results were consistent with those predicted in The Cancer Cell Line Encyclopedia (CCLE) database (Figure 7D and E). IHC sections from the Human Protein Atlas (HPA) database source showed that DDX56 is overexpressed at the protein level in cancerous tissues (Figure 7F), and we confirmed this result through our own samples (P < 0.001, Figure 7G and H).

Figure 7 Expression analysis of DDX56 in LIHC. (A) Differential expression of DDX56 between C1 and C2 clusters based on TCGA-LIHC databases (n = 374, ****P < 0.0001). (B) Expression level of DDX56 between the tumor and paired normal tissues in the TCGA-LIHC databases (n = 50, ***P < 0.001). (C) Differential expression of DDX56 between normal tissue and HCC tissue by using public datasets from Gene Expression Omnibus GSE36376 (n = 433, ****P < 0.0001) and GSE102079 (n = 257, ****P < 0.0001). (D and E) DDX56 protein expression in different HCC cell lines (**P < 0.01, ***P < 0.001). (F) Immunohistochemical analysis of DDX56 in HCC tissues and normal tissues based on Human Protein Atlas database (HPA). (G and H) Representative images of different immunohistochemical staining intensities for DDX56 based on our own HCC samples and statistical comparison of DDX56 expression levels (H-SCORE) in HCC tissues and normal tissues (n = 50, ***P < 0.001). (I) UMAP plots showing expression of DDX56 clusters. (J) The cell type assignment of each cluster is based on the DDX56 expression. (K) Correlation of the DDX56 expression levels with MKI67 based on TCGA-LIHC (n = 374). (L and M) Correlation between DDX56 and Ki67 based on immunohistochemical H-score calculation (n = 50). (N) Association of DDX56 mRNA expression with T/N/M stages, pathological stage, histologic grade and vascular invasion in HCC patients based on TCGA-LIHC (n = 374, NS: not significant, *P < 0.05, **P < 0.01, ***P < 0.001). (O) Analysis of differentially immune infiltration level with DDX56 high and low expression in TCGA-LIHC (NS: not significant, *P < 0.05, **P < 0.01, ***P < 0.001).

Single-cell analysis revealed that DDX56 are ubiquitously expressed by all cell types (Figure 7I and J). Considering the strong correlation between helicase phenotype and cell proliferation, we further found that the level of DDX56 was highly positively correlated with the level of MKI67 based on TCGA-LIHC calculation (R = 0.395, P < 0.001). Since both DDX56 and Ki67 are nuclear proteins, we found that they also converge in expression at the protein level by sequential comparison of IHC staining results (Figure 7L and M, R = 0.426, P = 0.002). Subsequently, we analyzed the relationship between DDX56 and clinicopathological parameters and found that the level of DDX56 increased with tumor grade as well as with vascular invasion (Figure 7N, P < 0.01). Finally, our calculations showed that the level of DDX56 significantly and negatively correlated with the abundance of most immune cells, including TGD cells, natural killer (NK) cells, dendritic cells (DC) and cytotoxic cells, among others (Figure 7O, P < 0.01). This part of the results indicates that DDX56 is highly expressed in hepatocellular carcinoma tissues and is a proliferative marker that is associated with immunosuppression.

Discussion

Helicases are essential for all living organisms, and their main function is to unpackage the genes of the organism.³ They act as linear molecular engines, specifically, moving directionally along the nucleic acid phosphodiester backbone and using the energy generated by ATP hydrolysis to separate two annealed nucleic acid strands, including DNA, ribonucleic acid (RNA) or RNA-DNA hybrids.^38–40 Thus, whenever cellular processes involving the separation of nucleic acid strands, such as DNA replication, transcription, translation, recombination, DNA repair and ribosome biogenesis, are closely associated with helicases, ultimately determining cell fate.^41,42

We examined the genetic and transcriptional profiles of 112 helicase molecules and found extensive mutational covariation among them, and not only that we also observed that the vast majority of helicases are overexpressed in almost all common cancer types. One of the essential features of tumor cells is their rapid proliferation accompanied by DNA damage response and chromosomal instability, therefore enhanced activity of resolvases is currently considered necessary for tumor cells to cope with endogenous or exogenous stress in order to maintain cellular homeostasis.^43–45 In addition, considering the frequent DNA damage caused by chemotherapy or radiotherapy, it has also been shown that DNA helicases play a key role in resistance to conventional clinical anti-cancer therapies, and Nicola J Curtin’s work suggested that the activity level of individual DNA damage response genes, including the decapping gene, can be used as a predictive biomarker for personalized cancer therapy.^46,47

To further understand the significance of helicases in LIHC, we pooled the expression profiles of helicase-encoding genes based on TCGA-LIHC and classified LIHC patients into two clusters (C1, C2) by consensus clustering, and patients in C1 versus C2 showed large differences in survival and clinical characteristics. Notably, levels of almost all helicases were observed to be significantly higher in C1 than in C2, and mutational profiling showed that TP53, a “star” gene currently found to prevent tumor cell division, induce apoptosis, and repair normal DNA damage,^48,49 was mutated much more frequently in C1 than in C2. Further, we found far more mutational co-occurrence events in C1 than in C2. This together with the above results suggested that C1 patients may be characterized by stronger genomic instability compared with C2. GSEA based on differential genes of two subtypes showed that KRAS signaling and hedgehog pathway activity were significantly upregulated in C1. To effectively explore deeper biological mechanisms, we identified key co-expressed gene modules in differential genes by MEGENA algorithm, and enrichment analysis showed that co-expressed gene modules are closely associated with cell cycle and division. Subsequently, we further calculated to obtain 10 hub genes, which have potential promotive effects on cell cycle. Interestingly, the levels of these 10 hub genes, including CDCA5, KIF2C, KIFC1, KIF4A, HJURP, TOP2A, BUB1B, KIF18B, KIF23, and GINS1, were significantly higher in tumor tissues than in paired normal tissues, and a similar phenomenon was observed between C1 and C2. These results show their potential oncogene function and possible relevance to nuclear division. We note that most of them belong to the kinesin family, and kinesins are now widely recognized as ATPases that move along microtubule filaments involved in the transport of cellular cargo and provide support for several cellular functions, including mitosis and meiosis.^50–52 The biological significance of kinesins in tumors is currently receiving increasing attention, and it has been demonstrated that KIF2C is a direct target of the Wnt/β-catenin pathway and mediates crosstalk between Wnt/β-catenin and mTORC1 signaling by interacting with TBC1D7.⁵³ In addition to these, the value of CDCA5 as a cell cycle-associated protein for nuclear division needs no elaboration,^54,55 while other hub genes including TOP2A, HJURP, BUB1B, and GINS1 have been widely reported to be risk factors for poor cancer prognosis, and their potential mechanisms are widely involved in interfering with DNA topology, mitotic checkpoint signaling, spindle assembly, and other biological processes that maintain genomic and chromosomal stability.^56–58

To distinguish in more detail the biological differences between the two helicase subtypes, we proceeded to examine the expression characteristics of these 10 hub genes at the single cell level. Considering the significant difference in the degree of mitotic activity between the two subtypes, we used the tumor proliferation index Ki67 as a reference to analyze the relationship between 10 hub genes and cell division. We observed a cluster of hepatocellular carcinoma cells at the G2M stage, enriched in Ki67-encoding genes, indicating that they were in a proliferative phase. Interestingly, most of the single-cell localizations of hub genes are consistent with Ki67, especially TOP2A, a cell proliferation marker that exerts a strong anti-apoptotic effect by frequently amplifying in tumor tissues.⁵⁹ In addition, KIFC1 was also observed to be abundantly expressed on proliferating cells, and some investigators observed reduced levels of Ki67 and PCNA after knockdown of KIFC1 in Hep3B and SNU-475 cell lines, tentatively confirming a potential positive correlation between KIFC1 and HCC proliferation.⁶⁰

The above work confirms the existence of two subtypes with significantly different proliferative activity in LIHC patients, and since this subtype is established based on the helicase expression profile, it is clear that the helicase activity itself may be a potentially critical factor. We identified three independent prognostic factors, DHX58, DDX43, and DDX56, through the development of the LASSO-COX model. Survival analysis showed that DDX56 was strongly associated with poor prognosis and that DDX56 was overexpressed in most tumor types.^61,62 We then focused on the significance of DDX56 at LIHC and further confirmed the high expression of DDX56 in 2 independent datasets and obtained consistent results at the protein level. Several studies have indicated that DDX56 interacts with the mitotic regulator serine/threonine phosphoprotein phosphatase (PP) 1 and the oncoprotein DEK to regulate tumor progression, and its high expression has been reported in breast and colorectal cancers.^62–64 Based on this, we compared the correlation between DDX56 and MKI67 and a significant positive correlation was observed at both transcriptome and protein levels. Finally, we also found that DDX56 levels negatively correlated with the abundance of most immune cells, which may be due to immunosuppressive effects caused by malignant proliferation of tumor cells.

Conclusions

In conclusion, our study highlights the relationship between helicase subtypes and proliferative features of LIHC. In conclusion, our study highlights the relationship between helicase isoforms and proliferative features of LIHC. These observations may be beneficial for conventional chemotherapy-based treatment of LIHC patients, especially for apoptosis-inducing and cell cycle-blocking-related interventions. We also identified DDX56 as an independent prognostic factor in patients with LIHC, and DDX56 was overexpressed in HCC tissue and significantly positively correlated with Ki67, indicating a high proliferation rate of tumor cells. The main limitation of this study was the small sample of clinical cases. The number of clinical cases was limited, mainly because of the death cases during the follow-up. Besides this, changes in living systems involve multi-scale regulatory processes in time and space. Various intracellular response mechanisms are finely regulated through complex signaling pathways, and the different pathways are intertwined with each other. DDX56, which has not been previously reported in hepatocellular carcinoma. The underlying mechanisms merit further explorations.