Abstract

   The liver exhibits extensive circadian regulation among organs.
   Epidemiological studies have substantiated that disruptions in
   circadian rhythm constitute a risk factor for the oncogenesis of liver
   cancer. Nonetheless, the molecular underpinnings of how circadian
   dysregulation influences liver cancer progression remain elusive. Our
   research aims to elucidate these mechanisms and develop a predictive
   model for prognosis and treatment responsiveness. Our multi-omics
   analysis revealed extensive dysregulation of liver circadian genes
   (LCGs) in liver cancer. Employing machine learning algorithms, we
   pinpointed four pivotal dysregulated LCGs. Through the integration of
   single-cell, bulk, and spatial transcriptomics, we further elucidated
   the interconnections between LCGs dysregulation and the tumor
   microenvironment. In vivo and in vitro experiments demonstrated that
   RBM17, identified as a crucial dysregulated LCG, promotes the
   progression of liver cancer and cisplatin resistance by facilitating
   cancer stem cell phenotype. The circadian prognosis scores (CPS), based
   on these four genes, effectively reflected the prognosis of liver
   cancer patients and their responses to various therapeutic
   interventions. Mechanism of Action (MOA) analysis suggested that high
   CPS level may sensitize tumors to cell cycle-targeted therapies.
   Collectively, our findings provide new insights into the interplay
   between liver circadian gene regulation and liver cancer progression,
   and propose novel therapeutic targets for liver cancer.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12935-025-03843-6.

   Keywords: Biomarkers, RBM17, Liver cancer, Circadian genes

Introduction

   Liver cancer remains a significant global health challenge, with an
   estimated incidence of over 1 million cases by 2025. The most common
   form of liver cancer is hepatocellular carcinoma (HCC) and it accounts
   for approximately 90% of cases [[46]1]. Surgery is the primary
   treatment for liver cancer patients. However, for those with advanced
   or metastatic liver cancer, treatment options are very limited. The
   first-line therapy has predominantly included sorafenib and
   transarterial chemoembolization (TACE), which have yielded suboptimal
   median overall survival (OS) rates. Despite immune checkpoint
   inhibitors exhibiting considerable therapeutic promise, the response
   rates are generally modest [[47]2, [48]3]. Consequently, there is an
   urgent need to develop robust models predictive of the prognosis and
   therapy response in liver cancer. Additionally, exploring the
   underlying mechanisms will identify potential targets and provide a
   theoretical basis for drug development and clinical decision making.

   Humans exhibit robust circadian rhythms that regulate a wide range of
   physiological processes. The circadian clock is an evolutionarily
   ancient and highly conserved timing mechanism that generates ~ 24-hour
   rhythms at multiple levels, from cells to organs to whole organisms,
   and synchronizes behaviors [[49]4]. This intrinsic clock system enables
   living organisms to anticipate and respond to daily environmental
   changes by temporally coordinating gene expression profiles.
   Epidemiological studies suggest the potential involvement of the
   deregulated circadian clock in the pathophysiology of various
   malignancies. Of note, the liver is a highly metabolically active organ
   that is strongly influenced by feeding and fasting cycles, and
   circadian misalignment [[50]4, [51]5]. Disrupted circadian rhythms are
   widely recognized to be involved in the pathological processes of
   various liver diseases, including nonalcoholic fatty liver disease
   (NAFLD), non-alcoholic steatohepatitis (NASH), and liver fibrosis,
   which can ultimately contribute to the development of liver cancer
   [[52]6, [53]7].

   Predictive models based on circadian genes are noteworthy due to their
   importance in liver cancer. However, most of such models were developed
   upon core clock genes [[54]8–[55]10]. A subset of core clock genes,
   such as CLOCK and BMAL1, among others, form the interlocked
   transcription-translational feedback loops that drive the rhythmic
   expression of circadian output genes. It’s estimated that hundreds and
   even thousands of gene expression in mice is associated with the
   circadian rhythm [[56]4, [57]5, [58]11]. Moreover, multiple levels of
   genetic, epigenetic, and transcriptional regulation contribute to the
   molecular clockwork. Therefore, current core clock biomarkers may not
   characterize tumor subtypes well due to the molecular heterogeneity of
   liver cancer. Integrating multi-omics analyses, including
   transcriptomic, genomics, and methylation sequencing data, with
   machine-learning algorithms can help predict the comprehensive
   association between circadian genes and liver cancer prognosis and
   therapy.

   Despite their clear association, the exact mechanism underlying the
   tumorigenic effects of the circadian genes remains largely unknown.
   Circadian rhythm disorders may trigger or at least be involved in
   metabolic dysregulation, cell cycle control, angiogenesis, genomic
   instability, cancer stem cell (CSC) maintenance, immune evasion, drug
   sensitivity, etc [[59]12–[60]14]. However, many of these proposed
   mechanisms still lack experimental evidence and require further
   functional analyses and validation. Moreover, the variation of the
   circadian clock among different cell types and its intricate interplay
   with the immune system and microenvironment niche makes the mechanism
   research more challenging. Bulk technologies could merely reflect the
   merged signals from heterogeneous cell groups. Therefore, single-cell
   analyses are essential to further dissect activities of circadian
   signaling pathways in a cell-type manner and predict their effects on
   immunotherapy response and immune features in liver cancer.

   Here, we unveiled the extensive dysregulation of circadian genes and
   their intricate interactions with malignant pathways in liver cancer
   through a comprehensive multi-omics landscape analysis. We analyzed the
   dysregulation of pivotal circadian genes within the tumor
   microenvironment using single-cell datasets and machine learning
   algorithms. Additionally, our research explored the significant role of
   the key circadian prognostic gene, RBM17, in the regulation of CSCs,
   and confirmed its contributory effect on cisplatin resistance of liver
   cancer. We introduced novel circadian prognosis scores (CPS) that
   accurately reflect the prognosis and treatment response of liver
   cancer. Moreover, we proposed the potential therapeutic vulnerability
   of tumors with high circadian disruption to cell cycle-targeted drugs.

Materials and methods

Public data sources acquisition and preprocessing

   The bulk RNA expression profiles containing the corresponding clinical
   survival information of liver cancer were collected from The Cancer
   Genome Atlas (TCGA,
   [61]https://www.cancer.gov/ccg/research/genome-sequencing/tcga) and
   International Cancer Genome Consortium (ICGC,
   [62]https://dcc.icgc.org/), including TCGA-LIHC (training set) and
   ICGC-LIRI-JP (validation set). Raw counts were converted to transcripts
   per kilobase million (TPM), and log2 transformed (TPM + 1). The
   Sorafenib, TACE and anti PD1 treatment bulk RNA sequencing cohorts were
   download from the Gene Expression Omnibus (GEO,
   [63]https://www.ncbi.nlm.nih.gov/geo/) through accessing id
   [64]GSE109211, [65]GSE104580 and [66]GSE202069. Single-cell
   RNA-sequencing (scRNA-seq) data of [67]GSE149614 were obtained from GEO
   [[68]15].

Multi-omics analysis of liver circadian genes

   Liver circadian genes (LCGs) were inferred using the cyclic ordering by
   periodic structure (CYCLOPS) algorithm developed by Ron C. Anafi et al.
   [[69]16], and they have been compiled into the CircaDB
   ([70]http://circadb.hogeneschlab.org/human) by Pizarro, A et al.
   [[71]17]. The multi-omics profiles data, including transcriptomic,
   genomics, and methylation arrays, are all sourced from the paired
   TCGA-LIHC cohort. The methylation data of TCGA-LIHC cohort are
   downloaded from XNEA ([72]https://xena.ucsc.edu/). Similar to the
   previous study, the methylation probe that shows the strongest negative
   correlation with mRNA expression was identified as the representative
   gene methylation probe used in this study [[73]18]. The “maftools” R
   package was applied to analyze single nucleotide variants (SNVs) of
   circadian genes. Copy number variation analysis was conducted using
   GISTIC2, with copy number alterations defined as greater than or equal
   to 2 or less than or equal to -2. The “ActivePathway” R package
   [[74]19] was utilized to analyze multi-omics pathway enrichment of
   LCGs. In brief, we calculated the multi-omics differences between tumor
   samples and paired normal tissues, and integrated pathway enrichment
   analysis was performed using the “ActivePathways” function. The
   merge_method parameter was set to “fisher”, and the cutoff parameter
   was set to “0.1”. Pathway enrichment was based on the KEGG database and
   visualized using Cytoscape [[75]20].

Single-cell RNA-sequencing analysis

   The downloaded data was processed using R package “Seurat” [[76]21].
   After quality control, the data were integrated by “Harmony” method. A
   total of 30,248 cells (16,093 tumor tissue sourced and 14,155 normal
   tissue sourced) were analyzed, with an average of 7,735 transcripts
   detected per cell. We performed dimensionality reduction using Uniform
   Manifold Approximation and Projection (UMAP). Finally, the cells were
   annotated by classical markers. The copy number variations and stemness
   of tumor cells were predicted using R package “inferCNV” and
   “CytoTRACE” [[77]22]. The R package “Monocle” was used to calculate and
   visualize the pseudotime developmental trajectories of tumor cells
   [[78]23].The EMT signature and Proliferation signature were obtained
   from the MsigDB database ([79]https://www.gsea-msigdb.org) and
   calculated by the “AUCell” method. For drug prediction analysis, the
   “Beyondcell” R package were utilized [[80]24].

Spatial transcriptomics data processing

   For the ST-seq data analysis, we utilized the “Seurat” R package
   (v4.1.1). The data, obtained from
   HRA000437([81]https://ngdc.cncb.ac.cn) [[82]25], underwent
   “SCTransform” normalization. Additionally, the expression pattern of
   RBM17 was visualized using the “SpatialFeaturePlot” function in
   “Seurat”.

Selecting the key feature of liver circadian genes by machine-learning

   To reveal the potential prognosis and molecular mechanism of liver
   circadian genes (LCGs), we identified the optimal prognostic biomarkers
   among the LCGs using the LASSO Cox regression model with the “glmnet” R
   package in the TCGA-LIHC cohort. To select the most suitable lambda
   value and avoid overfitting, we applied ten-fold cross-validation. A
   total of 10 LCGs with nonzero coefficients were selected through
   10-fold cross-validation. Random Forest is a versatile and accessible
   machine-learning algorithm often used for classification and regression
   tasks. As it has no restrictions on variable conditions, it has higher
   accuracy, sensitivity, and specificity than decision trees in
   predicting continuous variables and obtains predictions without
   significant bias [[83]26]. Here, we used the R Package Random Forest to
   screen for the most significant LCGs associated with prognosis. Default
   parameters were used as they consistently produced stable results:
   ntree = 500. A total of the top 10 LCGs were selected based on the rank
   of “MeandecreaseAccuracy”. Finally, by combining Random Forest and
   Lasso regression, we took the intersection of the 2 subsets of feature
   genes selected from the above methods to ultimately determine the key
   feature genes.

Construction of circadian prognosis scores (CPS)

   The circadian prognosis scores (CPS) were constructed using regression
   coefficients based on IL18RAP, GNL2, RBM17, ZDHHC18. The CPS was
   formulated as follows:
   graphic file with name d33e411.gif

   Patients were stratified into high and low CPS groups based on the
   median CPS. The Kaplan-Meier method was employed to assess differences
   in survival between the high- and low-CPS groups. Additionally,
   time-dependent receiver operating characteristic (ROC) curve analysis
   was conducted to evaluate the accuracy of the CPS in predicting
   TCGA-LIHC prognosis. The robustness of the CPS was further validated in
   the ICGC cohort.

Pathway-enrichment and immune infiltration analysis

   For pathway enrichment, gene sets were obtained from the MsigDB
   database ([84]https://www.gsea-msigdb.org). The “GSVA” R package was
   applied to calculate pathway enrichment characteristics for single
   samples. The “IOBR” R package serves as an integrated tool for TME
   analysis [[85]27]. The immune infiltration was computed using the
   “deconvo_tme” function of the “IOBR”. Specific parameters were set to
   their default values without any modifications. Additionally, the MHC
   signature was obtained from the “IOBR”. These signatures were employed
   in estimating the relative scores for each sample using the “GSVA”
   method.

Clinical samples

   A total of 112 HCC tissue samples were collected with informed consent
   from patients who underwent hepatectomy at the Sun Yat-Sen University
   Cancer Center (Guangzhou, China). Complete clinic pathological and
   follow-up data were available, and patient demographic information is
   provided in Supplementary Table [86]S1. The use of clinical specimens
   in this study was approved by the Committee for Ethical Review of
   Research Involving Human Subjects at the Sun Yat-Sen University Cancer
   Center.

Immunohistochemistry (IHC) staining

   Immunohistochemistry (IHC) analysis was performed using a primary
   antibody against RBM17 (1:200 dilution, Cat#13918-1-AP, Proteintech).
   The IHC staining followed standardized protocols as described elsewhere
   [[87]28]. Each sample was scored as 1 (negative), 2 (weak positive), or
   3 (strong positive) based on the relative staining intensity. Samples
   exhibiting strong positivity were classified as the high expression
   group, while negative and weakly positive samples were categorized into
   the low expression group.

Immunofluorescence staining

   The cell slides were fixed with 4% paraformaldehyde, then were
   pre-incubated with 10% normal goat serum at room temperature for 30 min
   to reduce nonspecific reaction. Subsequently, the slides were incubated
   with rabbit monoclonal anti-CD133 (1:200 dilution), for 12 h at 4 °C.
   The slides were incubated with a secondary antibody (fluor 488 labeled)
   for 1 h at room temperature and stained with DAPI for 5 min. The
   proportion of CD133 + cells was calculated by counting the number of
   CD133-positive cells in randomly selected regions of the slide. The
   total number of cells was determined by counting DAPI-stained nuclei in
   the same regions. The ratio of CD133 + cells was calculated as the
   number of CD133-positive cells divided by the total number of cells in
   each field.

Cell culture

   The HCC cell line PLC-8024 (CRL-8024) and HepG-2 were acquired from the
   Institute of Virology of the Chinese Academy of Sciences (Beijing,
   China). 293T cells were obtained from National Collection of
   Authenticated Cell Cultures, Shanghai, China. All cell lines were
   authenticated by short tandem repeat profiling and routinely tested for
   the absence of mycoplasma contamination. Specifically, we used a
   PCR-based mycoplasma detection kit (GMyc-PCR Mycoplasma Test Kit,
   Cat#40601ES10, Yeasen) to test the cell culture supernatant at regular
   intervals. The cells were cultured in Dulbecco’s modified Eagle’s
   medium (DMEM) (Gibco), supplemented with 10% fetal bovine serum (Gibco)
   and 1% penicillin/streptomycin mixture (Gibco). All cell lines utilized
   in this investigation were incubated at 37 °C in a humidified incubator
   maintained at 5% CO[2].

Plasmids, lentivirus production and cell infection

   The RBM17-shRNA expression vectors and the scrambled shRNA control were
   constructed with the pLL3.7 plasmid (Addgene). The RBM17-shRNA target
   sequence is GATGAAGCAGTACGGATATTT. Virus was made by co-transfecting
   the pLL3.7 constructs along with helper plasmids from the Directional
   TOPO Expression Kit (Invitrogen)into 293T cells using HieffTrans
   Liposomal Transfection Reagent (Yeasen Biotech) according to the
   manufacturer’s recommendations. After transfection for 72 h, the viral
   supernatant was collected and used to infect the cell line.
   Subsequently, stable knockdown cell lines were obtained by further
   selection with puromycin (Gibco).

Cell viability assay and foci formation

   Cells were seeded into 96-well plates at a density of 2000 cells per
   well in 100 µl of growth medium. Subsequently, 10 µl of CCK8 solution
   (Cat# HY-K0301, MCE) was added to each well, followed by incubation at
   37 °C for 3 h. The optical density was measured at a wavelength of
   452 nm using a microplate reader over a period of 5 days. Triplicate
   repeats were performed to evaluate variance and statistical
   significance. The data points were represented as the average of
   triplicate wells. As for foci formation, cells were plated in 6 well
   plate. 7–12 days later cells were washed with PBS three times, fixed
   with 75% ethanol for 30 min, and dyed purple by gentian violet. The
   colony formation rates were calculated by dividing the number of foci
   by the total number of cells initially seeded, expressed as a
   percentage.

Transwell and wound healing assay

   For the migration assay, 2 ~ 5 × 10^4 cells were suspended in 100 µL
   DMEM and plated onto the upper chamber (Corning) inserts. Meanwhile,
   700 µL DMEM supplemented with 10% FBS was added to the lower chamber in
   a 24-well plate. After 48-hour incubation at 37 °C, the migrated cells
   were fixed and dyed. For the invasion assay, chambers were uniformly
   covered with 60 µL Matrigel (Cat# 354234, BD Biosciences) diluted with
   DMEM (1:8). After 2-hour incubation at 37 °C, cell suspension was
   plated onto the upper chamber, with 700 µL DMEM supplemented with 10%
   FBS adding to the lower chamber. Finally, the invasive cells were dyed.
   For the wound healing assay,2 ~ 5 × 10^4 cells were plated onto 6-well
   plates. When the confluence was up to 90%, the cell monolayer was
   scratched with a 1 mL pipette tip. After removing detached cells, the
   remaining cells were cultivated in medium without FBS. Cellular
   migration was quantified using Image J.

Sphere formation assay

   Cells were plated as single cells into 96-well plates at 100 cells/well
   with 2 mL DMEM/F12 (Hyclone) containing 1% B27 supplement
   (Thermofisher), 1% N2 supplement (Thermofisher), 20 ng/mL EGF
   (Thermofisher) and 10 ng/mL bFGF (Thermofisher). After 7 days of
   growth, spheres were counted. Sphere formation rate = sphere-formed
   cells/seeded cells.

Western blot

   Samples were lysed with Pierce RIPA buffer (Thermo Fisher Scientific,
   USA) and quantified using Bio-Rad Protein Assay (Bio-Rad). After being
   mixed with loading dye (Bio-Rad), the protein was denatured (100℃,
   15 min). Denatured proteins were separated by SDS-PAGE and transferred
   onto a PVDF membrane (Millipore, USA). After blocking for 1 h, the
   membranes were incubated with primary antibodies (RBM17, Cat#
   13918-1-AP, Proteintech,; GAPDH, Cat# 2118, CST; CDK2, Cat# 2546, CST;
   CyclinE2, Ca# 4132, CST; CD133, Cat# 18470-1-AP, Proteintech) and then
   incubated with peroxidase-conjugated secondary antibodies.
   Immunoreactive bands were visualized using ECL (Bio-Rad) and exposed to
   autoradiograph film.

Cell cycle and apoptosis flow cytometry

   Cells were grown to a density of 70% in 6-well dishes and collected for
   staining with propidium iodide (PI, Cat# A601112-0100, Sangon)
   according to the manufacturer’s instructions. For the apoptosis assay,
   the Annexin V-Alexa Fluor 488/7-AAD Cell Apoptosis Detection Kit (Cat#
   40313ES50, Yeasen Biotech) was utilized, and the procedures were
   conducted according to the manufacturer’s instructions. A CytoFLEX
   cytometer (Beckman Coulter) was used to measure changes in cell cycle.
   FlowJo was utilized to analyze the generated flow cytometry data.

Animal experiment

   The Institutional Animal Care and Use Committee at Southern University
   of Science and Technology has reviewed and approved all animal
   experiments. Four- to five-week-old male BALB/c-nu/nu mice were reared
   under specific pathogen-free conditions at the Laboratory Animal Centre
   of Southern University of Science and Technology. Afterwards, the mice
   were randomly divided into four groups (nc, nc + Cisplatin, shRBM17,
   shRBM17 + Cisplatin). In the drug resistance assay, PLC-8024 cells
   stably transfected with lenti-sh-RBM17 or lenti-sh-NC were injected
   subcutaneously into BALB/c-nu/nu mice (4 × 106 cells, 5 mice per
   group). When the tumor volume reached approximately 50 mm^3. Mice were
   treated with cisplatin (3 mg/kg every three days, intraperitoneally),
   while control mice were treated with saline. Tumor size was measured
   every 3 days. Tumor volume was estimated using the formula: V = length
   × width^2/2. All mice were housed under a controlled 12-hour
   light/12-hour dark cycle. Euthanasia was performed during the standard
   daytime working period.

Statistical analysis

   Data processing, analysis, and visualization were carried out using R
   packages (version 4.30; [88]https://www.bioconductor.org/) or Graphpad
   Prism software (version 8.0.1; [89]https://www.graphpad.com/). To
   compare differences between two groups, we employed Student’s t-test or
   Mann–Whitney U test. Correlation analysis between variables was
   performed using Pearson or Spearman test. Survival analysis was
   conducted utilizing the Kaplan–Meier approach and log-rank test, with
   the survival and “Survminer” packages. Uni- and multivariable Cox
   regression models were constructed to compute hazard ratios (HRs) and
   identify independent prognostic factors. A significance level of
   P < 0.05 was considered statistically significant. Results are
   expressed as the mean ± standard error of the mean.

Results

Multi-omics analysis reveals dysregulation of liver circadian genes in liver
cancer

   Our overall research workflow is presented in (Fig. [90]1A). A total of
   378 Liver circadian genes (LCGs) were gathered from CircaDB, which are
   genes rhythmically expressed in the normal liver. Transcription
   profiles of 378 LCGs showed extensive disruption in liver tumors
   compared with normal livers from the TCGA-LIHC dataset (Fig. [91]2A).
   We integrated multi-omics data to investigate whether genetic,
   epigenetic, and transcriptional regulation contributes to the
   disordered expression of LCGs. Genetic mutations of LCGs are not
   frequent (less than 5%) in liver cancer, as illustrated in
   (Fig. [92]2B). Regarding copy number variation (CNV) characteristics, a
   considerable portion of LCGs displayed prevalent copy number
   amplifications (Fig. [93]2C). Methylation characteristics of LCGs were
   also investigated. The top 15 genes with upregulated methylation and
   the top 15 genes with downregulated methylation in tumors were
   visualized (Fig. [94]2D). Multi-omics correlation analysis reveals a
   complex interplay in the dysregulation of LCGs. The overall methylation
   associated with LCGs expression shows a negative correlation, while a
   positive correlation is observed with copy number variations. BMAL1,
   CLOCK, NR1D1 and NR1D2 are considered as core transcription factors in
   circadian regulation [[95]29–[96]32]. Their expression levels
   positively correlate with the overall LCGs expression, with CLOCK
   demonstrating the strongest correlation among them (Fig. [97]2E).
   Multi-omics integration pathway enrichment analysis indicates that
   dysregulation of LCGs is associated with pathways of cancer, cell
   proliferation, metabolism, endothelium-interstitial transition, innate
   immunity, and adaptive immunity. LCGs mRNA dysregulation plays a
   predominant role in most pro-malignant pathways (Fig. [98]2F and
   Supplementary Table [99]S2).

Fig. 1.

   [100]Fig. 1
   [101]Open in a new tab

   Overall research workflow of this study. (A). The flowchart of the
   study. TF, Transcription factor; LCGs, Liver circadian genes; TME,
   Tumor microenvironment; CSC, Cancer stem cell; TCGA, The Cancer Genome
   Atlas; ICGC, International Cancer Genome Consortium

Fig. 2.

   [102]Fig. 2
   [103]Open in a new tab

   Multi-omics analysis reveals dysregulation of liver circadian genes in
   liver cancer. (A). Expression heatmap of 378 LCGs in TCGA-LIHC
   datasets. (B). Mutations of LCGs in the TCGA-LIHC cohort. Each column
   corresponds to an individual sample, and each row to a specific genes.
   Mutation type is indicated by a distinct color. (C). CNV variant
   frequencies of LCGs in TCGA-LIHC. The height of the column represents
   the frequency of change. (D). The dotplot of alterations in LCGs
   methylation. (E). The ridge plot depicting the multi-omics correlation
   of circadian genes. Blue indicates a negative correlation, red
   indicates a positive correlation, and the dashed line marks the median
   correlation. (F). The ActivePathway enrichment plot for dysregulated
   LCGs. Node size represents the number of genes enriched in the pathway,
   while node color corresponds to the respective omics dimension. CNV,
   Copy number variation; TF, Transcription factor

Identifying critical LCGs for liver cancer prognosis using machine learning

   Given the extensive crosstalk between LCGs disruption and pro-malignant
   pathways, we further employed machine learning techniques to identify
   LCGs associated with liver cancer prognosis. The Random Forest
   algorithm and LASSO regression are two machine-learning techniques used
   for identifying the most significant predictors in a dataset. Our
   approach combined these two algorithms to enhance the
   representativeness and stability of feature selection. Ultimately,
   IL18RAP, RBM17, GNL2, and ZDHHC18 were pinpointed as the most pivotal
   dysregulated circadian genes in predicting liver cancer prognosis
   (Fig. [104]3A). Kaplan-Meier curves reveal ZDHHC18, GNL2, and RBM17 as
   adverse prognostic factors in liver cancer, while IL18RAP emerges as a
   favorable prognostic factor in both TCGA and ICGC cohorts
   (Fig. [105]3C). Consistently, ZDHHC18, GNL2, and RBM17 exhibit
   significant upregulation in tumor tissue compared to the normal bulk
   gene expression level, while IL18RAP shows downregulation
   (Fig. [106]3B).

Fig. 3.

   [107]Fig. 3
   [108]Open in a new tab

   Identifying critical LCGs for liver cancer prognosis using machine
   learning. (A). The strategy for selecting key LCGs using machine
   learning. (B). Box plots of gene expression in TCGA adjacent normal and
   tumor tissues. (C). The responding Kaplan-Meier curve for ZDHHC18,
   IL18RAP, GNL2, and RBM17. The first four figures correspond to results
   obtained from the TCGA dataset, while the next four figures pertain to
   results from the ICGC cohort. The Kaplan–Meier curves were analyzed
   using the log-rank test. TCGA, The Cancer Genome Atlas; The Cancer
   Genome Atlas; ICGC, International Cancer Genome Consortium

Dissecting critical LCGs activities in tumor microenvironment using
single-cell RNA transcriptome

   To further dissect the activities of LCGs in various cellular
   compositions within the tumor microenvironment of liver cancer, we
   carried out scRNA-seq analysis. After implementing quality control
   measures, we integrated scRNA-seq data from ten tumor tissues and eight
   paired normal liver tissues obtained from the [109]GSE149614 dataset
   (Supplementary Figure [110]S1A). The cells were then annotated as “B
   cells” (B), “Fibroblasts” (Fibro), “Endothelial cells” (Endo),
   “Hepatocytes”, “Myeloid cells” (Myeloid), and “T/NK cells” (T/NK)
   (Fig. [111]4A-B and Supplementary Figure [112]S1B). In exploring the
   expression of four key LCGs within the tumor microenvironment, we noted
   a widespread upregulation of RBM17 and GNL2 in tumor-sourced
   hepatocytes. RBM17 was also found to be enriched in tumor-derived B
   cells and fibroblasts. However, no significant changes were observed in
   myeloid and endothelial cells from different tissue sources, suggesting
   that RBM17 may exert a potentially tumor-promoting effect in certain
   cell types (Fig. [113]4C). Despite the low expression levels of ZDHHC18
   and IL18RAP on a single-cell basis, likely due to dropout, our analysis
   identified specific expression pattern of IL18RAP in T/NK cells,
   especially those from normal tissues (Fig. [114]4D and Supplementary
   Figure [115]S1C). Given the pivotal role of T cells in the tumor
   microenvironment, we subclustered T/NK cells for further investigation
   (Supplementary Figure [116]S1D). This led to the clustering of T cells
   into several subsets: cytotoxic CD8 T cells, exhausted CD8 T cells, T
   regulatory cells (Tregs), naive T cells, CD4 T cells, natural killer
   (NK) cells, and a small fraction of other cells (Fig. [117]4E).
   Notably, a significant upregulation of IL18RAP was observed in
   cytotoxic CD8-GZMK cells and, to a lesser extent, in NK cells
   (Fig. [118]4F). This suggests a potential role of IL18RAP in enhancing
   immune surveillance. Supporting our hypothesis, our bulk RNA-Seq
   analysis showed a significant positive correlation between IL18RAP and
   immune-related pathways, as well as with CD8T/NK cells (Supplementary
   Figure [119]S1E).

Fig. 4.

   [120]Fig. 4
   [121]Open in a new tab

   Dissecting critical LCGs activities in tumor microenvironment using
   single-cell RNA transcriptome. (A). Dotplot of marker gene expression
   in different cell types. (B). UMAP graph of cell types and sample
   sites. (C). Expression violin-plot of core circadian prognostic genes,
   with green indicating adjacent normal tissues and red indicating tumor
   samples. (D). Expression Dotplot of IL18RAP in different cell types.
   (E). UMAP graph of different T cell subtypes. (F). Gene relative
   expression heatmap of different T cell sub-cluster. cyCD8T, cytoxic
   CD8T cell; ex CD8T, exhausted CD8T cell; Treg, regulatory T cell; NK,
   Natural killer cells; B, B cell; Endo, Endothelial cell

Multi-omics analysis and experimental validation unravel the pro-malignant
role of RBM17 in liver cancer

   The gene RBM17 was selected for further exploration due to its
   significantly increased expression in tumor-sourced hepatocytes. We
   initially confirmed the significant upregulation of RBM17 in
   leading-edge section, tumor section, and portal vein tumor thrombus
   section in publicly available spatial transcriptomic slices.
   (Supplementary Figure [122]S2A). Subsequent immunohistochemical (IHC)
   staining validated the significant association between RBM17 expression
   and worse overall survival in our in-house HCC cohort. (Fig. [123]5A).
   We then examined the expression of RBM17 across different tumor-sourced
   hepatocyte clusters using public scRNA-seq data [[124]15]. The
   hepatocytes were subgrouped into 10 clusters, with RBM17 showing the
   most significant upregulation in cluster 2 (Fig. [125]5B and
   Supplementary Figure [126]S2B). The inferCNV tool was utilized to
   calculate the copy number variation(CNV) scores for the various tumor
   clusters (Figure [127]S2C). The findings revealed that cluster 2
   exhibited a high CNV score, which is frequently linked to tumor
   malignancy. The epithelial-mesenchymal transition (EMT) score and
   proliferation capability were assessed using the AUCells method.
   Consistently, cluster 2 showed the highest scores among the 10 tumor
   clusters (Fig. [128]5C). The Monocle pseudotime analysis further
   revealed that cluster 2 resides at the trajectory endpoint, indicating
   a progressive increase in RBM17 expression along tumor progression
   (Fig. [129]5D). In the bulk RNA sequencing level, single sample gene
   set enrichment analysis (ssGSEA) was conducted in the TCGA-LIHC cohort,
   which revealed a significant positive correlation between RBM17
   expression and pathways regulating cell proliferation, stemness, and
   EMT (Supplementary Figure [130]S2D). In view of the above evidence, the
   necessity of RBM17 for the malignant functions of liver cancer cells
   was further examined. In vitro experiments demonstrated that
   shRNA-mediated silencing of RBM17 reduced cell viability and the
   frequency of foci formation in human HCC cell lines PLC-8024 and HepG2
   (Fig. [131]5E). In addition, wound healing assays and transwell
   experiments revealed that knocking down RBM17 significantly inhibited
   the (EMT) phenotype in PLC-8024 and HepG2 cells (Fig. [132]5F-H). In
   summary, the analyses and experiments consistently demonstrate the
   significant role of RBM17 in liver cancer progression.

Fig. 5.

   [133]Fig. 5
   [134]Open in a new tab

   Multi-omics analysis and experimental validation unravel the
   pro-malignant role of RBM17 in liver cancer. (A). Representative images
   of immunohistochemistry (IHC) staining showing high and low expression
   levels of RBM17 in tumor tissues of HCC (left panel). The Kaplan-Meier
   curve of different RBM17 expression groups (right panel). The
   Kaplan–Meier curves were analyzed using the log-rank test(n = 112).
   (B). UMAP graph of subsetted tumor cell clusters (Top panel). Dimplot
   of RBM17 expression level (Bottom panel). (C). Boxplots of the CNV
   scores, EMT scores, and proliferation scores across different clusters.
   (D). Monocle pseudo-temporal plot. The top panel represents the pseudo
   time. The middle panel illustrates the positions of different cell
   clusters within the pseudo trajectory, with the C2 cluster (high RBM17
   expression) located at the trajectory’s terminus. The bottom panel
   represents the relative expression levels of RBM17 in pseudo time. (E).
   Cell proliferation curve by cck8 assay (left panels, n = 3) and foci
   formation assay (right panels, n = 3). (F). Representative image of
   wound healing assay of the indicated cells. The corresponding barplot
   on the right illustrates the wound healing rate(%). (G). Representative
   image of transwell assay of PLC-8024(n = 3). (H). Representative image
   of transwell assay of HepG-2(n = 3).
   *P < 0.05,**P < 0.01,***P < 0.005,****P < 0.001, The data are presented
   as the means ± standard deviations (SD) of three independent
   experiments. The P value was determined using the student t-test

RBM17 regulates the cancer stem cell properties in liver cancer

   Previous studies have suggested an association between disrupted
   circadian rhythms and cancer stem cells (CSCs) [[135]33]. Our
   multi-omics analysis also revealed a correlation between RBM17 and
   liver cancer stemness and differentiation. Next, we aim to investigate
   whether RBM17 affects the stemness of liver cancer cells. CytoTRACE
   analysis showed a significant overlap between RBM17 expression and
   poorly differentiated cells, suggesting an association between RBM17
   and CSCs (Fig. [136]6A). Furthermore, analysis of CSC signature based
   on TCGA-LIHC cohort revealed that patients with high RBM17 expression
   had higher mRNA stemness index (mRNAsi) scores (P < 0.01), indicating
   increased stemness (Fig. [137]6B). Based on these results, we
   hypothesize that RBM17 may help in regulating the stemness of tumor
   cells. In support of this, spheroids formation assays demonstrated that
   silencing RBM17 significantly reduced the number of spheroids formed by
   CSCs (Fig. [138]6C). Pseudo-bulk analysis of RBM17 highly-expressing
   cells revealed the upregulation of many CSCs markers, including CD133
   (Fig. [139]6D). Consistently, the protein expression level and
   immunofluorescence signals of CD133 were significantly reduced after
   RBM17 depletion (Fig. [140]6E-F). CSCs exhibit high proliferative
   properties. Flow cytometry analysis revealed that the loss of RBM17 led
   to cell cycle arrest at the G1 phase, accompanied by reduced expression
   of G1 phase-promoting proteins CDK2 and Cyclin E2 (Fig. [141]6G-H).
   Overall, these results underscore the critical role of RBM17 in
   maintaining CSCs and consequently promoting the progression of liver
   cancer.

Fig. 6.

   [142]Fig. 6
   [143]Open in a new tab

   RBM17 regulates the cancer stem cell properties in liver cancer. (A).
   UMAP plot of CytoTRACE predicted cell order (more to less
   differentiated), indicated by color from blue to red (left panel). UMAP
   plot of RBM17 expression levels (right panel). (B). Boxplot of mRNA
   stemness index (mRNAsi). (C). Representative image of spheroids formed
   assay and quantification of the indicated cells. (D). Lollipop plot of
   CSCs marker fold change. (E). Inhibition of CD133 by RBM17 silencing
   was verified by western blotting of the indicated cells. GAPDH was used
   as the loading control. (F). Immunofluorescence analysis of CD133
   expression in PLC-8024 and HepG-2 cells. The nuclei were stained with
   DAPI (blue). (G). Flow cytometry detection of cell cycle using PI
   staining in PLC-8024 and HepG-2(n = 3). (H). RBM17 and G1/S-related
   proteins, including CyclinE2, CDK2 were detected with WB.
   *P < 0.05,**P < 0.01,***P < 0.005,****P < 0.001, The data are presented
   as the means ± standard deviations (SD) of three independent
   experiments. The P value was determined using the student t-test

RBM17 suppression enhances cisplatin sensitivity in liver cancer

   Drug resistance is a defining characteristic of cancer stem cells(CSCs)
   [[144]34]. Given RBM17’s essential role in maintaining cancer stemness,
   we explored its relationship with drug resistance. Using the
   “Beyondcell” algorithm, we predicted the heterogeneity in drug
   responses of liver cancer cells to cisplatin and sorafenib at the
   single-cell level. Our findings revealed that cells expressing high
   levels of RBM17 significantly overlapped with those exhibiting
   resistance to cisplatin (Fig. [145]7A-B). This was illustrated by a
   notable increase in the BSC values for cisplatin resistance in tumor
   cells with high RBM17 expression, in comparison to other tumor cells.
   However, no significant variance was observed in response to sorafenib
   between the two cell groups (a higher BSC value indicates greater drug
   resistance) (Fig. [146]7B). In addition, the “Oncopredict” tool was
   employed to predict the sensitivity of various patients to these drugs,
   based on data from the TCGA-LIHC cohort. Within this cohort, RBM17
   expression levels were positively correlated with the IC50 values of
   sorafenib (R = 0.28) and cisplatin (R = 0.52) (Fig. [147]7C),
   suggesting a potential involvement of RBM17 in liver cancer drug
   resistance. Furthermore, analyses based on the TACE cohort consistently
   demonstrated that RBM17 serves as a predictive marker for chemotherapy
   drug response (Fig. [148]7D). It is noteworthy that our results
   consistently demonstrate a significant association between RBM17
   expression and cisplatin resistance across different computational
   methods. In vitro assays showed that RBM17 depletion reduced the
   survival of PLC-8024 and HepG-2 cell lines across various cisplatin
   concentration gradients (Fig. [149]7E). Apoptosis analysis further
   demonstrated that knocking down RBM17 significantly enhanced
   cisplatin-induced apoptosis in PLC-8024 (P < 0.001) (Fig. [150]7F). To
   corroborate the role of RBM17 in drug resistance in vivo, we
   established a cell line-derived xenograft (CDX) in nude mice.
   Consistent with our in vitro findings, RBM17 knockdown markedly
   inhibited tumor growth and amplified the inhibitory effects following
   cisplatin treatment in vivo. (Fig. [151]7G).

Fig. 7.

   [152]Fig. 7
   [153]Open in a new tab

   RBM17 suppression enhances cisplatin sensitivity in liver cancer. (A).
   UMAP plot of tumor cells and the definite cell population (left). UMAP
   plot for BCS Score mapping. The red color reflects resistance to
   Cisplatin (right). (B). Boxplot of BSC value in different RBM17 status
   cell groups. (C). Scatter plot of the correlation between predicted
   IC50 and the expression of RBM17. (D). Boxplot of RBM17 expression in
   different TACE response groups (left). The Area Under Curve (AUC) of
   RBM17 reflects the TACE response. (E). The survival curves based on the
   CCK8 assay(n = 3). (F). Flow cytometry apoptosis analysis depicting the
   enhancement of cisplatin-induced apoptosis in PLC-8024 following RBM17
   silencing(n = 3). (G). Photos of excised subcutaneous Cell line-derived
   xenograft tumors and quantification of tumor weight and volume.
   *P < 0.05,**P < 0.01,***P < 0.005,****P < 0.001.The data are presented
   as the means ± standard deviations (SD) of three independent
   experiments. The P value was determined using the student t-test

Development of circadian prognosis scores to enhance liver cancer prognosis
prediction

   We explored the substantial correlation between the dysregulation of
   four circadian genes, notably RBM17, and the prognosis of liver cancer.
   To quantitatively evaluate the prognosis of liver cancers, we devised
   the circadian prognosis score (CPS). The CPS distribution and survival
   status across in TCGA and ICGC cohorts were delineated in the top panel
   (Fig. [154]8A). For the TCGA cohort, the average AUC values for 1-year,
   3-year, and 5-year prognosis predictions were 0.77, 0.72, and 0.70,
   respectively (Fig. [155]8A). Correspondingly, for the ICGC cohort, the
   average AUC values for 1-year, 2-year, and 3-year survival were found
   to be 0.73, 0.71, and 0.70, respectively (Fig. [156]8A). To augment the
   practical applicability of the CPS, we amalgamated independent
   prognostic factors and formulated a nomogram. This nomogram serves as a
   quantitative instrument, enabling clinicians to ascertain the mortality
   risk associated with liver cancers (Fig. [157]8B). The calculation of a
   total score for each patient is achieved by summing the scores
   attributed to each prognostic parameter. It is noted that a higher
   total score in indicative of a worse prognosis for the patient.
   Kaplan-Meier curves elucidated the superior survival discrimination
   capability across distinct NomoRisk categories (Fig. [158]8C).
   Additionally, a calibration plot affirmed that the performance of the
   nomogram closely parallels that of an ideal model (Fig. [159]8D),
   further validating its efficacy and reliability in the clinical
   prognosis of liver cancer.

Fig. 8.

   [160]Fig. 8
   [161]Open in a new tab

   Development of circadian prognosis scores to enhance liver cancer
   prognosis prediction. (A). Survival status, time distribution and model
   genes expression of patients in the high circadian prognosis scores
   (CPS) and low CPS groups (Top panels). Time-dependent receiver
   operating characteristic (ROC) curve analysis showing the area under
   curve values (AUC) for overall survival of TCGA and ICGC cohorts
   (Middle panels). Kaplan–Meier survival curves for liver cancer patients
   according to different risk types in TCGA and ICGC cohorts (Bottom
   panels). (B). The comprehensive nomogram for predicting probabilities
   of liver cancer patients with 1-, 3- and 5-year OS in the TCGA dataset.
   (C). Kaplan–Meier survival curves for liver cancer patients according
   to different nomoRisk types in TCGA cohort. (D). The calibration plots
   for predicting liver cancer patients with 1-, 3- and 5-year OS in the
   TCGA cohort. Nomogram-predicted probability of survival is plotted on
   the x-axis; actual survival is plotted on the y-axis. The Kaplan–Meier
   curves were analyzed using the log-rank test. CPS, Circadian prognostic
   scores; AUC, Area under the ROC Curve; OS, Overall survival; TCGA, The
   Cancer Genome Atlas; The Cancer Genome Atlas; ICGC, International
   Cancer Genome Consortium

Different genomic mutation landscapes and tumor-associated pathways among
high- and low-risk groups

   We further explored the molecular and cellular characteristics of the
   two different risk groups defined by CPS. It is well-established that
   genomic mutations are pivotal in the initiation and progression of
   tumorigenesis [[162]35]. Analysis of the single-nucleotide variation
   (SNV) landscape revealed distinct gene mutation profiles of different
   CPS populations (Fig. [163]9A). Notably, the mutation frequencies of
   TP53, along with SP1 and RESF1, were found to be higher in patients
   with high CPS, while TOGARAM2 showed a higher mutation frequency in
   patients with low CPS (P < 0.01) (Fig. [164]9B). Specifically, the TP53
   gene, a critical tumor suppressor gene, showed mutations in 79 out of
   182 patients with high CPS, while only 20 out of 179 patients with low
   CPS exhibited mutations (P < 0.005) (Fig. [165]9B). Mutation site
   analysis uncovered extensive alterations in the TP53 gene body region
   among individuals with high CPS, including the hotspot mutation R249S
   (Fig. [166]9C). Additionally, co-mutation analysis revealed a
   widespread co-occurrence of mutations among individuals within the high
   CPS group (Fig. [167]9D). Similarly, copy number variation analysis
   indicated that the high CPS group exhibited more copy number
   amplifications (CNA) and loss of heterozygosity (LOH). Pathway
   enrichment analysis identified distinct pathway enrichment features in
   individuals with high CPS, implicating DNA replication, cell cycle
   regulation, DNA damage repair, and metabolism (Fig. [168]9F). In
   conclusion, these results highlight significant disparities in gene
   mutations and pathway activation among different CPS populations,
   offering valuable insights into the potential roles of circadian gene
   dysregulation in liver tumorigenesis.

Fig. 9.

   [169]Fig. 9
   [170]Open in a new tab

   Different genomic mutation landscapes and tumor-associated pathways
   among high- and low-risk groups. (A). Mutations of top 20 genes in the
   TCGA-LICH cohort. Each column represents each sample, and each row
   represents each gene. Mutation type is marked by a distinct color. The
   left panel represents the high CPS groups and the right panel
   represents the low CPS groups. (B). Forest plots of mutation bias in
   high CPS and low CPS groups. (C). Lollipop plots of mutation sites in
   different CPS groups. (D). Exclusive/co-occurrence mutation event
   heatmap. (E). Copy number amplification (CNA) fraction and Loss of
   heterozygosity (LOH) fraction altered boxplot in different CPS groups.
   (F). KEGG pathway enrichment heatmap in different CPS groups.
   *P < 0.05,**P < 0.01,***P < 0.005,****P < 0.001, The P value was
   determined using the student t-test. CPS, circadian prognostic scores;
   TMB, tumor mutational burden; CNA, copy number amplification; LOH, loss
   of heterozygosity

Distinct immune microenvironments and responses to therapy among high- and
low-risk groups

   An increasing body of research highlights a significant interplay
   between circadian rhythms and the tumor immune microenvironment (TME).
   To investigate the association between CPS and TME, we employed a
   comprehensive approach, integrating algorithms such as EPIC, MPC, IPS,
   xCell, and Estimate, to generate detailed immune infiltration profiles
   for different CPS subgroups (Fig. [171]10A). The results revealed
   pronounced immune cell infiltration in the low CPS subgroup, in
   contrast to the observed immune exhaustion within the high CPS
   subgroup. Furthermore, the MHC signature, as determined through GSVA,
   demonstrated elevated levels in the low CPS populations, indicative of
   enhanced immune functionality (Fig. [172]10B and C). These findings
   suggest that a low CPS is associated with a “hot” immune subtype,
   characterized by with robust immune activity, which may be more
   amenable to immunotherapy. Supporting this hypothesis, out analysis
   utilizing the TIDE algorithm identified patients likely to respond to
   immunotherapy, who predominantly exhibited lower CPS values. A notable
   positive correlation was established between CPS and TIDE score, with a
   higher TIDE score suggesting a decreased likelihood of benefiting from
   immunotherapy (R = 0.42) (Fig. [173]10D-F). Employing CPS to
   differentiate immunotherapy response yielded an AUC value of 0.65
   (Fig. [174]10G). Consistency in these findings was further validated in
   a cohort of liver cancer patients treated with anti-PD-1 therapy
   (Supplementary Figure [175]S3A). These results indicate that CPS
   defines distinct immune subtypes and may be used to guide immunotherapy
   decisions. Despite the promising potential of immunotherapy for
   treating liver cancer, sorafenib and TACE therapy remain the mainstream
   treatment modalities. Analysis of patient cohorts undergoing sorafenib
   treatment revealed a significant increase in CPS among non-responders,
   with an AUC value of 0.88 for sorafenib response prediction. Similarly,
   within cohorts receiving TACE therapy, non-responders exhibited
   elevated CPS levels, with an AUC value of 0.70 for response prediction
   (Fig. [176]10H-I). To identify potential therapeutic agents for high
   CPS patients, the GDSC2 datasets and “Oncopredict” algorithm were
   employed to ascertain the correlation between CPS and sensitivity to a
   panel of 192 drugs. Subsequent ranking identified 10 drugs with a
   correlation coefficient (“Rho”) less than − 0.4, indicating potential
   efficacy. These agents include ML323_1629, Sepantronium, Tozasertib,
   BMS.345541_1249, MK.1775_1179, Daporinad, Ulixertinib, Paclitaxel,
   GDC0810_1925, and Wee1.Inhibitor_1046 (Fig. [177]10J). An in-depth
   analysis of the mechanism of action (MOA) of these identified compounds
   revealed a predominant targeting of cell cycle-related mechanisms,
   underscoring a putative vulnerability of high CPS patients to
   interventions targeting the cell cycle (Fig. [178]10K). Consistently,
   GSEA analysis revealed significant activation of cell cycle-related
   pathways in patients with high CPS (Supplementary Figure [179]S3B).

Fig. 10.

   [180]Fig. 10
   [181]Open in a new tab

   Distinct immune microenvironments and responses to therapy among high-
   and low-risk groups. (A). Immune infiltration landscape of different
   CPS populations. The barplot on the left represents the correlation
   between CPS and immune cells. (B). Boxplot of MHCI scores calculated by
   GSVA method. (C). Boxplot of MHCII scores calculated by GSVA method.
   (D). Boxplot of CPS values in different ICB treatment response
   populations. (E). Scatterplot of correlation between Tide score and
   CPS. (F). Sankey diagram of CPS and immunotherapy response. (G). The
   Area Under Curve (AUC) of CPS reflects ICB response (H). Boxplot of CPS
   in different Sorafenib treatment response populations (left panel). The
   Area Under Curve (AUC) of CPS reflects the Sorafenib response (right
   panel). (I). Boxplot of CPS in different TACE treatment response
   populations (left panel). The Area Under Curve (AUC) of CPS reflects
   the TACE response(left panel). (J). Correlations rank plot of CPS and
   Estimated drugs IC50. Blue represents a negative correlation, while red
   represents a positive correlation. The dashed position represents the
   selected threshold. (K). MOA Dotplot for Screening Drugs Negatively
   Correlated with CPS. Pearson correlation analysis was performed to
   assess the correlations. Differences between the two groups were
   analyzed using a student t-test. CPS, Circadian prognostic scores; MHC,
   major histocompatibility complex; NR, non-response; R, response; AUC,
   Area under the ROC Curve; IC50, half-maximal inhibitory concentration

Discussion

   Our current study unveiled the multi-omics dysregulation of circadian
   rhythm in liver cancer, highlighting the intricate crosstalk between
   rhythm transcription factors, promoter methylation, copy number
   variations, and circadian dysregulation. Utilizing the ActivePathway
   approach, we constructed a multi-omics functional landscape of
   disordered circadian genes. The findings reveal extensive modulation of
   pathways implicated in cancer, EMT, cell proliferation, and
   immune-related pathways, as a consequence of circadian genes
   disruption. These insights underscore the pivotal role of circadian
   genes disturbance in the advancement of liver cancer.

   Research employing humanized mouse models has elucidated that circadian
   disruption directly reprograms the transcriptome landscape of NASH and
   HCC, thereby accelerating the transition from HCC onset to the
   progression of HCC [[182]7]. Despite these advances, pinpointing the
   specific critical circadian genes responsible for the advancement of
   liver cancer remains a considerable challenge.

   In our study, RBM17, GNL2, ZDHHC18, and IL18RAP were identified as key
   dysregulated circadian genes through a machine learning approach.
   Furthermore, single-cell sequencing analysis revealed the dysregulation
   of these key circadian genes within the tumor microenvironment.
   Notably, our analysis highlighted the distinct expression of IL18RAP in
   cytotoxic CD8^+ T cells, alongside its subsequent downregulation in the
   hepatic tumor milieu. This observation may elucidate the mechanism
   behind immune suppression induced by circadian dysregulation, offering
   insights into the intricate relationship between circadian gene
   expression and immune function within the tumor microenvironment.
   Another critical circadian signature gene, RBM17, known as a
   splicing-related RNA-binding protein, was further validated in our
   in-house HCC samples. RBM17 expression is markedly elevated in HCC
   relative to normal liver tissue and correlates with adverse prognostic
   outcomes. Our comprehensive multi-omics investigation, integrating,
   bulk, and single-cell transcriptomic datasets, indicates that the
   upregulation of RBM17 is associated with increased CNV scores, enhanced
   cell proliferation, EMT, poor differentiation, and augmented stemness
   characteristics. Experimental studies, both in vitro and in vivo,
   elucidate that RBM17 significantly contributes to the proliferative and
   migratory properties of liver cancer cells. Furthermore, RBM17 is
   implicated in promoting CSC stemness and CD133 expression, alongside
   facilitating CSC-associated traits, especially in cisplatin resistance.
   Our study expands on the understanding of RBM17’s role in mediating
   cisplatin resistance in HCC. Notably, CD133 has been well-established
   as a critical regulator of cisplatin resistance in cancers, with CD133⁺
   cells demonstrating marked resistance to cisplatin treatment through
   mechanisms such as downregulation of PTEN and activation of the
   PI3K/AKT/mTOR signaling pathway [[183]36, [184]37]. RBM17-mediated
   regulation of CD133 expression may therefore explain its involvement in
   modulating cisplatin sensitivity. The investigation into the role of
   RBM17 across various cancers, though limited, reveals a consistent
   pattern. Specifically, RBM17 is implicated in controlling cell
   proliferation in HCC, glioma, and hypopharyngeal carcinoma
   [[185]38–[186]40]. Moreover, it enhances CSC properties in colorectal
   cancer and acute myeloid leukemia, and influences the radiosensitivity
   and cisplatin sensitivity of lung and hypopharyngeal cancer cells,
   respectively [[187]38, [188]41–[189]43]. Additionally, disruptions in
   the circadian rhythm have been observed to significantly exacerbate
   cancer cell dissemination and lung metastasis in breast cancer, further
   enhancing tumor cell stemness and the tumor-initiating potential, while
   inducing an immune-suppressive metastatic tumor microenvironment
   [[190]44]. Collectively, our findings furnish both bioinformatic and
   experimental substantiation linking RBM17 with liver tumorigenesis and
   elucidate a mechanistic connection between disrupted circadian rhythms
   and the enhanced stemness of cancer cells.

   The Circadian Prognosis Scores (CPS) were developed to reflect
   circadian dysregulation in liver cancers. This approach may
   significantly optimize risk stratification, highlighting the role of
   liver circadian genes as robust prognostic classification factors.
   Furthermore, CPS is capable of precisely predicting the therapeutic
   responses of liver cancers to mainstream treatments such as sorafenib
   and TACE, in addition to immunotherapy. Moreover, CPS delineates
   distinct immune subtypes and pathway-enriched subtypes. Despite the
   resistance to TACE, sorafenib, and immunotherapy in cases presenting
   high CPS, analysis of the mechanism of action (MOA) reveals the
   therapeutic promise of drugs targeting the cell cycle for patients with
   elevated CPS. Supporting this, evidence derived from study on
   regenerating mouse livers suggests that the circadian clock regulates
   the expression of cell cycle-related genes, which in turn modulate the
   expression of the active cyclin B1-Cdc2 kinase, a key regulator of
   mitosis [[191]45]. These findings underscore the vital interplay
   between circadian rhythm and cell cycle modulation, suggesting a
   pronounced vulnerability of tumors with high CPS to cell cycle-targeted
   therapeutics.

   The evolutionary history of liver cancer is rooted in its genome
   [[192]46]. Among the genes, TP53 is the top mutated gene in HCC, with
   mutation frequencies ranging from 18 to 35.2% (average 25.9%)
   [[193]47]. Intriguingly, our study indicates a notable bias towards
   high circadian disruption associated with TP53 mutations. This suggests
   an unexplored role of TP53 in circadian disruption in liver cancer.
   Supporting our findings, p53^−/− mice have a shorter period length that
   lacks stability, along with impaired photo-entrainment to a light
   pulses in a free-running condition [[194]48]. Therefore, the intricate
   relationship between TP53 mutations and circadian disruption merits
   further investigation.

   This study is subject to several limitations. First, the LCGs analyzed
   were identified through computational methods using the CYCLOPS
   algorithm, which infers rhythmicity without requiring time-of-day
   annotations. The adoption of this approach was necessitated by the
   formidable practical and ethical challenges associated with acquiring
   time-resolved human liver samples. Consequently, these LCGs are
   circadian-associated rather than rhythmically oscillating genes within
   the human sample datasets. Further experimental validation and
   phase-specific profiling are required. Secondly, although the
   prognostic significance of CPS was corroborated across multiple
   independent cohorts, our study lacks an internal cohort for validation
   purposes. Additionally, the evaluation of CPS’s utility in
   immunotherapy was constrained by limited data, relying on the TIDE
   algorithm and a small cohort of HCC patients treated with PD-L1
   inhibitors, which underscores the need for validation in a broader
   patient cohort. Lastly, despite suggesting that high CPS populations
   may benefit from cell cycle-targeted therapies, this hypothesis lacks
   validation through organoid models derived from patients.

   In summary, our findings lay a foundational understanding of the
   interplay between liver circadian genes and liver cancer progression,
   paving a new avenue for identifying relevant therapeutic targets.

Conclusions

   We constructed a CPS that is capable of predicting the prognosis and
   response to various treatments, including immune therapy. Furthermore,
   bioinformatics and functional evidence confirmed that RBM17, within the
   CPS, may play crucial regulatory roles in liver cancer. This positions
   RBM17 as a potential therapeutic target and prognostic biomarker for
   liver cancer.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [195]Supplementary Material 1^ (95.5KB, xlsx)
   [196]Supplementary Material 2^ (15.9KB, docx)
   [197]Supplementary Material 3^ (95.5KB, xlsx)
   [198]Supplementary Material 4^ (1.2MB, docx)

Acknowledgements