Identification of a Novel Circadian Rhythm-Related Signature for Predicting Prognosis and Therapies in Hepatocellular Carcinoma Based on Bulk and Single-Cell RNA Sequencing

,


Introduction
Hepatocellular carcinoma (HCC) stands as one of the most prevalent malignant liver cancers, accounting for 75%-85% of cases and ranking as the fourth leading cause of cancerrelated deaths globally [1].Hepatitis B virus (HBV) or hepatitis C virus (HCV) infection, afatoxin-contaminated food, high alcohol consumption, smoking, obesity, nonalcoholic fatty live disease (NAFLD), and type 2 diabetes are the main risk factors that contribute to HCC development [2,3].Over the past decade, advancements in HCC management have signifcantly improved the overall survival and quality of life for patients [4].Commonly, patients with early-stage HCC are often considered for surgical resection or transplantation.Tose in the intermediate stage typically undergo the standard of care, involving transarterial chemoembolization.For patients in the advanced stage, systemic therapies with frst-line and second-line drugs, such as atezolizumab (anti-PD-L1 antibody), bevacizumab (anti-VEGF antibody), sorafenib, and lenvatinib (tyrosine kinase inhibitors, TKIs), are accepted treatment options [5].However, there are limited clinical benefts because of the frequent development of resistance [6].Te tumor microenvironment (TME) of HCC is a complex and structured mixture characterized by abnormal angiogenesis, chronic infammation, and dysregulated extracellular matrix (ECM) remodeling, which contribute to the hypoxia, immunosuppressive, and acidic microenvironment [7,8].Te complex TME mediates the aggressive tumorigenesis, metastasis, therapeutic resistance, immune evasion, and recurrence of HCC [8][9][10].Terefore, exploring the molecular mechanisms in the regulation of TME of HCC and then identifying novel therapeutic strategies that improve survival are the next major challenges in HCC.
Te circadian clock is a complex cellular mechanism that sustains self-perpetuating oscillations with a 24-hour periodicity to control several cyclic physiological processes, and disrupting circadian rhythm results in numerous physiological disorders and diseases, including cancer [11,12].Emerging evidence has demonstrated that circadian clockcontrol metabolism is a hallmark of cancer, and circadian rhythm can be found as the novel therapeutic target [13][14][15][16].Circadian rhythm plays a pivotal role in tumorigenesis and tumor progression by orchestrating various biological processes within cancer cells.Tese processes include cell proliferation, apoptosis, DNA repair, and metabolism.In addition, circadian rhythm infuences the TME by shaping the cellular properties and interactions among various components, such as tumor-associated macrophages (TAMs), myeloid-derived suppressor cells (MDSCs), neutrophils, dendritic cells, T cells, natural killer (NK) cells, cancer-associated fbroblasts (CAFs), and endothelial cells [17,18].Given the tissue-specifc properties of circadian clocks, the liver serves as a physiological hub for circadian regulation [19,20].Te dysfunction of the hepatic circadian clock is implicated in the regulation of various liver functions, encompassing synthetic and metabolic processes related to glucose, lipids, bile acids, amino acids, and more [21].Te circadian clock genes contribute to tumor growth and aggressive [20,22,23].Alterations of the circadian clock genes are associated with the prognosis of HCC patients [24].Te abovementioned evidence has demonstrated that circadian rhythm plays an important role in hepatocarcinogenesis.
In the present study, we integrated the single-cell RNAseq (scRNA-seq) data and bulk RNA-seq data to comprehensively identify the prominent cell types and elucidate the cell-to-cell communication in TME that was infuenced by circadian rhythm.Moreover, selecting the prognostic circadian rhythm-related genes (CRRGs) based on scRNA-seq data and bulk RNA-seq data analyses to construct a risk model for survival prediction and the biological functions, immune infltration characteristics, and therapeutic responses associated with CR-risk score were investigated.Our fnding elucidated the TME characteristics by circadian rhythm infuences and provided novel potential therapeutic options for HCC treatment.

Screening the Diferentially Expressed Circadian Rhythm-Related Genes (CRRGs)
. "Limma" package in R [25] was performed to screen the diferentially expressed genes (DEGs) between HCCs and nontumors with the thresholds of |log 2 (fold change, FC)| > 0.585 and P value of <0.05.Te CRRGs further were identifed by overlapping the DEGs and 85 CRRGs.

2.3.
Calculating the CR Score.CR score for each sample in the TCGA-LIHC cohort was calculated using single-sample gene set enrichment analysis (ssGSEA) through the "GSVA" package in R [26].Te "survminer" package in R [27] was performed to select the optimal cutof value according to the CR score of each sample, the survival time, and the survival status.Te patients were divided into high-and low-CR score groups according to the optimal cutof value.Te overall survival analysis was performed between high-and low-CR score groups.

Construction of a Weighted Gene Coexpression Network Analysis (WGCNA).
Te "WGCNA" package in R was used to construct a coexpression network and identify the CR score-related modules and genes [28].Te samples in the TCGA-LIHC cohort were normalized, and the outlier samples were removed.Te soft threshold power (β) was selected to ensure the network was scale-free, and the genes in the frst quartile of variance were calculated by a power function.Te adjacency matrix was then transformed into a topological overlap matrix (TOM), and the corresponding dissimilarity (1-TOM) also was calculated.Te dynamic tree-cut method was used to identify the module by hierarchically clustering genes.A deepSplit value of 2 and a minimum size cutof of 30 were selected as the distance measure for constructing the dendrogram.Te signifcant modules were screened based on the correlation between CR score and modules by Pearson's correlation test.Te gene expression profles with module eigengenes (Mes) were 2 European Journal of Cancer Care defned as a module membership (MM), and the correlation between outer features and gene expression profles was defned as the gene signifcance (GS).Te hub genes were identifed that are located in the modules with the highest MM and highest GS values.
2.5.Single-Cell RNA-Seq Data Analysis.Te normalized scRNA-seq data of the GSE156625 cohort was obtained from the GEO database, which comprised a total of 57, 25414 HCCs and 1 normal sample [29]."Scanpy," a scalable Python-based package [30], was used for downstream analysis.Unifed manifold approximation and projection (UMAP) was used for dimensionality reduction and cell clustering visualization [31].Te "AUCell" package in R [32] was performed to calculate CR scores for each cell type, and the cells were divided into high-CR score and low-CR score cell populations based on the median value of AUCell value.Te diferential gene expression between high-CR and low-CR score groups was identifed using "scanpy.tl.rank_ge-nes_groups" [30].Te "GSEApy" package was performed for biological function analysis [33].

Construction and Validation of a Circadian Rhythm-Related Signature for Predicting Prognosis. Te intersected
CRRGs between the DEGs from scRNA-seq data analysis and the hub genes from WGCNA were screened and visualized by a diagram analysis.A univariate Cox analysis was used to select the prognostic CRRGs in HCC.Ten, the least absolute shrinkage and selection operator (LASSO) regression was performed using the "glmnet" package [34] to shrink the list of genes and to obtain the risk coefcients strongly linked to prognosis, and a risk model was further constructed.Te risk score was calculated as follows: risk score = sum (each gene expression x the corresponding coefcient).HCC patients in TCGA-LIHC, LIRI-JP, and GSE14520 cohorts were divided into high-risk and low-risk groups with the median value of risk score.Te overall survival curves were drawn using the Kaplan-Meier method by the "survminer" R package [35], and the performance of the risk model was evaluated by the area under the curve (AUC) values of the receiver operating characteristic (ROC) curves which were generated using "survivalROC" package in R [35].In addition, GO and KEGG pathway enrichment analyses between high-risk and low-risk groups were performed using the "GSEApy" package in R [30].Reverse transcription was performed with Prime-Script RT reagent Kit (DBI® Bioscience, Shanghai, China) according to the manufacturer's instructions.For real-time PCR analysis, the resultant cDNA products were amplifed using SYBR Green qPCR Master Mix (DBI, Shanghai, China) in triplicates.Primer sequences of the genes analyzed were GAPDH: 5′-TGTTCGTCATGGGTGTGAAC-3′ and 5′-ATGGCA TGGACTGTGGTCAT-3′; RPL29 : 5′-ACACCACACACA ACCAGTCC-3′ and 5′-GCATTGTTGGCCTGCATCTT-3′; PFKFB3: 5′-GATGCCCTTCAGGAAAGCCT-3′ and 5′-GAACACTTTTGTGGGGACGC-3′; RPS7: 5′-CCAAGC GAAATTGTGGGCAA-3′ and 5′-CCTTGCCCGTGAGCT TCTTA-3′; SLC6A6: 5′-CCCAGGCTCTCTGAAATGGG-3′ and 5′-AGGAGCATGGCGAATGGAAA-3′; and RPLP2: 5′-CGACCGGCTCAACAAGGTTA-3′ and 5′-GGCTTT ATTTGCAGGGGAGC-3′.GAPDH was used for normalization of the expression levels of each gene.Te relative expression was quantifed by using the 2 −ΔΔCt method.

Construction of a Predictive Nomogram.
Te risk score and clinical characteristics (such as age, sex, and clinical grade) were subjected to multivariate Cox regression to obtain the independent risk factor, and the relevant forest plot was drawn using the "ezcox" package in R [36].A predictive nomogram was constructed using the "regplot" package in R [37] based on the independent risk factors.ROC curves were used to evaluate the performance of the predictive model.A calibration plot was drawn to calculate the efciency of the predictive model using the "rms" package [38].A decision curve was drawn to assess the ability of the predictive model by the "dcurves" package [39].

Estimation of the Immune Cell Infltration in the High-CR Risk and Low-CR Risk Groups.
A total of 28 immune cellrelevant gene sets were constructed according to the previous study [40].ssGSEA was performed using the "GSVA" package in R [26] to calculate the immune infltration score for each sample in the TCGA-LIHC cohort, and the differences between the high-risk and low-risk groups were tested using the Kruskal-Wallis test.Moreover, the diferential expression of immune checkpoint-related genes between high-risk and low-risk groups was evaluated using the Kruskal-Wallis test.

Evaluation of the Response to Immunotherapy and Chemotherapy in HCC.
Tumor immune dysfunction and exclusion (TIDE) was used to predict immune checkpoint blockade (ICB) response [41].Immunophenoscore (IPS) also can be used to predict the response to immunotherapies including CTLA-4 and PD-1 blockers [40].Here, the "tidepy" function in the Python package was performed to European Journal of Cancer Care evaluate the TIDE score of each sample in the TCGA-LIHC cohort [41], and the "IOBR" R package was performed to evaluate the immunophenoscore (IPS) of each sample in the TCGA-LIHC cohort [42].Te signifcant diferences in TIDE score and IPS between high-risk and low-risk groups were detected using the Kruskal-Wallis test.Moreover, the half-maximal inhibitory concentration (IC 50 ) for patients with HCC based on the Genomics of Drug Sensitivity in Cancer (GDSC2) database (https://www.cancerrxgene.org) was calculated by the "oncoPredict" package in R, which was performed to predict drug response [43].Te diferences in log 2 (IC 50 ) between high-risk and low-risk groups were tested using the Kruskal-Wallis test.
2.12.Statistical Analysis.Statistical analyses were performed for all experiments with the GraphPad Prism software (Version 8.0, San Diego, CA).Results are presented as mean ± SD from at least 3 independent experiments.Te statistical diferences were calculated by using the Student's t-test.* P < 0.05 versus the control group, * * P < 0.01 versus the control group, and * * * P < 0.001 versus the control group.S2).Ten, a total of 40 diferentially expressed CRRGs were obtained (Figures 1(c) and 1(d) and Table S3).Circadian rhythm plays an important role in the regulation of complex physiological activities in HCC [44].Terefore, the CR score for each sample in the TCGA-LIHC cohort was calculated using ssGSEA, and patients were distributed into high-and low-CR score groups based on the optimal cutof value.Te patients were divided into high-and low-CR score groups according to the optimal cutof value, and the results suggested that patients in the high-CR score showed poorer survival time than those in the low-CR score group (Figure 1(e)).Tese fndings suggested that dysregulation of circadian rhythm was associated with the diferential prognosis in HCC.

Identifcation of the Key Modules and Hub Genes Related
to Circadian Rhythm.WGCNA was constructed to identify the key modules and genes that were related to the circadian rhythm in HCC [45].After normalization and removing the outlier samples (Figure 2(a)), the coexpression network was constructed with the β-value as 14, and the scale-free R 2 was equal to 0. Te normalized scRNA-seq data (GSE156625) were performed using UMAP dimensionality reduction to identify 28 Louvain clusters using the "Scanpy" package (Figure 3(a)).Ten, 28 Louvain clusters were annotated into 11 major cell types based on the data source article, including hepatocytes, fbroblasts, bipotent cells, endothelial cells, B cells, myeloid, CD4+ cells, CD8+ cells, natural killer (NK) cells, Tregs, and mast cells (Figure 3(a)).Te GO annotation enrichment analysis based on the DEGs of each population indicated that distinct cell clusters were involved in diferent biological processes (Figure 3(b)).For example, B cells were associated with the regulation of humoral immune response, complement activation, and immune efector process.CD4+ T cells, CD8+ T cells, and NK cells were involved in posttranslational protein modifcation.CD4+ T cells and Tregs were involved in receptor-mediated endocytosis.CD8+ T cells, NK cells, and Tregs were involved in platelet degranulation.NK cells and Tregs were involved in regulated exocytosis.Ten, to explore the CRRG expression characteristics in HCC, the CR score of each cell was calculated using the AUCell R package (Figure 3(c)).All cells were distributed into the high-CR score and low-CR score groups based on the median value of AUCell value; hepatocyte, bipotent cells, fbroblasts, and endothelial cells were distributed into the low-CR score group; and most immune cells were distributed into the high-CR score group (Figure 3(c)).A total of 2792 DEGs (2649 upregulated and 143 downregulated) were identifed between the high-CR score and low-CR score groups (Table S5).GO and KEGG pathway enrichment analyses based on the DEGs between the high-CR score and low-CR score groups are exhibited in Figures 3(d) and 3(e), and upregulated DEGs are involved in cytoplasmic translation and cotranslational protein targeting the membrane.Te downregulated DEGs were associated with prostanoid and prostaglandin metabolic processes.Tese fndings suggested that immune cells in the TME were more likely to be infuenced by the circadian rhythm.

Construction and Validation of a Circadian Rhythm-
Related Signature.To explore the prognostic values of CRRGs in HCCs, a total of 44 CRRGs were obtained by overlapping CRRGs from WGCNA and scRNA-seq analyses (Figure 4(a) and Table S6).Univariate Cox analysis indicated that 21 CRRGs were associated with prognosis (Figure 4(b)).Ten, fve CRRGs (RPL29, PFKFB3, RPS7, SLC6A6, and RPLP2) were selected as the prognostic signature using the LASSO regression analysis (Figures 4(c)-4(e)).Ten, the CR-related risk scores were calculated based on the prognostic gene expression and corresponding coefcients, and HCC patients were classed into high-risk and low-risk groups according to the median value of risk scores both in training cohort (TCGA-LIHC, Figures 4(f

qPCR Validation of the Expression of Circadian Rhythm-
Related Signature.To validate the expression of the circadian rhythm-related signature (RPL29, PFKFB3, RPS7, SLC6A6, and RPLP2) in HCC, qPCR was conducted to examine the mRNA expression of RPL29, PFKFB3, RPS7, SLC6A6, and RPLP2.As shown in Figures 5(a found that the PFKFB3, RPS7, RPL29, RPLP2, and SLC6A6 were upregulated in HCC samples compared with controls.

Development of a Predictive
Nomogram.We incorporated the clinical characteristics into the univariate and multivariate Cox models to identify the independent risk factors in HCC, and the results suggested that risk score was selected as the independent risk factor of HCC (Figures 6(a  6 European Journal of Cancer Care 3.6.Biological Function Analyses.We also investigated the biological function between high-risk and low-risk groups using GSEA.Te GO analysis indicated that CR risk was involved in the cellular protein metabolic process, focal adhesion, steroid hydroxylase activity, and cell-substrate junction (Figures 7(a)-7(c)).Te KEGG pathway enrichment analysis indicated that CR risk was linked to several metabolism-related pathways, such as fatty acid, tryptophan, tyrosine, glycine, serine, threonine, retinol, and bile acid metabolisms, and other enrichment pathways including RNA transport, DNA repair, and Wnt-beta catenin signaling pathway (Figure 7(d)).Furthermore, the hallmark pathway enrichment analysis revealed that CR risk was connected to several tumorigenesis pathways, such as the Wnt-beta catenin signaling pathway, p53 pathway, IL-2/STAT5 signaling pathway, mTORC1 signaling pathway, and PI3K/AKT/ mTOR signaling pathway (Figure 7(e)).

Description of CR-Risk Score Relevant Tumor Immune Infiltration Landscape
We also described the immune cell landscape between highrisk and low-risk groups using ssGSEA.As shown in Figures 8(a) and 8(b), we found that most adaptive immune cells, such as the activated dendritic cells, CD56dim NK cells, immature dendritic cells, MDSC, macrophage, mast cell, NK T cells, and plasmacytoid dendritic cells, were signifcantly enriched in the high-risk group than in the low-risk group.
Moreover, the innate immune cells also revealed a similar outcome, and the activated B cells, activated CD4+ T cells, central memory CD4+ T cells, central memory CD8+ T cells, efector memory CD4+ T cells, immature B cells, memory B cells, regulatory T cells, T follicular helper cells, Type17 T helper cells, and Type2 T helper cells also increased in the high-risk group than in the low-risk group (Figures 8(c) and 8(d)).Tese results were consistent with the scRNA-seq data analysis, in which CR played an important role in the regulation of the tumor environment.

Prediction of the Response to Immunotherapy and Chemotherapy in High-Risk and Low-Risk
Groups.Next, we explored the immune checkpoint-related genes (ICRGs) expression between high-risk and low-risk groups, results showed that the upregulated ICRGs (BTLA, BTN2A1, CD160, CD209, CD226, CD27, CD276, CD28, CD40LG, CD47, CD70, CD80, CD86, CD96, CTLA4, HAVCR2, HLA-DRB1, ICOS, IDO1, LAG3, LGALS9, PDCD1, PDCD1LG2, TIGIT, TNFRSF14, TNFRSF18, TNFRSF4, TNFRSF9, TNFSF18, TNFSF4, and TNFSF9) were found in high-risk group compared with low-risk groups (Figures 9(a) and 9(b)).Ten, we investigated the TIDE value and IPS, and the results indicated that higher TIDE value and lower IPS were observed in the high-risk group than in the low-risk group (Figures 9(c) and 9(d)).Tese results implied that high-risk scores were associated with a negative response to immunotherapy.We also investigated the distinction of European Journal of Cancer Care chemotherapy sensitivity between the high-risk group and the low-risk group, as shown in Figures 9(e) and 9(f ), lower IC 50 values of MK-1775 (adavosertib), paclitaxel, pevonedistat, ULK1_4989, vinblastine, and vinorelbine in the highrisk group than the low-risk group.Te results indicated that more chemotherapy sensitivity was found in high-risk scores compared with the low-risk score group.

Discussion
Te circadian clock, an endogenous timekeeper system, consists of both master and peripheral clock genes, orchestrating various biological processes.Disruption of the circadian rhythm detrimentally impacts physiology and poses global health threats, contributing to proliferative, metabolic, and immune diseases [46,47].Circadian rhythmregulated metabolism is a novel hallmark cancer and involves hepatocarcinogenesis, progression, metastasis, treatment outcomes, recurrences, and survival [48][49][50][51].
Increasing research suggests that circadian rhythm disruption is associated with the risk and prognosis of colorectal cancer [52].Moreover, circadian rhythm disruption accelerates cancer growth, worse survival, and chemoresistance in pancreatic cancer [53].Most studies conducted        on the infuences of circadian rhythm so far have focused on hub circadian clock genes that observed expression changes in the experimental detection results.However, a large number of the CRRGs that contribute to tumorigenesis and development are still unknown.Tus, in the present study, we comprehensively and intensively investigated the hub CRRGs and their regulatory mechanisms in HCC.First, we identifed the diferentially expressed CRRGs in HCC and demonstrated that the diferential expression of CRRGs is involved in survival.Terefore, we further screened the 252 hub genes that are associated with circadian rhythm.Meanwhile, we also investigated the landscape of the TME based on scRNA-seq data.We found that 11 major cell types, including hepatocytes, fbroblasts, bipotent cells, endothelial cells, B cells, myeloid, CD4+ cells, CD8+ cells, NK cells, and Tregs, were identifed, and those genes were enriched in diferent biological processes.With the CR score calculation, all cells were clustered into high-CR score and low-CR score groups according to the median value of the CR score, and hepatocyte, bi-potent cells, fbroblasts, and endothelial cells were distributed into the low-CR score group, and most of the immune cells (B cells, myeloid, CD4+ cells, CD8+ cells, NK cells, and Tregs) were distributed into the high-CR score group.Tese data indicated that immune cells were strongly infuenced by the circadian rhythm.Cancer-immunity cycle is a circulation system in TME that is composed of several major steps such as cancer cell antigen release and presentation, priming and activation of efector immunity cells, tracking and infltration of immunity to tumors, and elimination of cancer cells [54].It has been found that both the innate and adaptive arms of immunity in the TME possess circadian rhythmicity and have been linked to antitumor response [55][56][57][58][59][60][61][62].Our results supported the abovementioned research fndings that immune cells in the TME are regulated by circadian rhythmicity.
Te abovementioned fnding reminded us that circadian rhythm is involved in tumor growth, the outcome of treatment, and survival.Tus, we further investigated the prognostic values of CRRGs in HCC.By integrating scRNA-seq data and bulk RNA-seq data analyses, fve CRRGs (RPL29, PFKFB3, RPS7, SLC6A6, and RPLP2) were selected and validated as the prognostic signature in HCC.Moreover, the experimental results also demonstrated the signifcant diferential expression of fve prognostic CRRGs.HCC patients were divided into high-risk and lowrisk groups based on the median value of risk score.We also found a high-CR risk score linked to several metabolismrelated pathways and canonical cancer-related pathways, such as the Wnt-beta catenin signaling pathway, p53 pathway, and PI3K/AKT/mTOR signaling pathway.Besides, a high-CR risk score not only promoted the increased proportions of innate and adaptive immune cells but also negatively associated with immunotherapeutic responses and positively associated with some chemotherapeutic drugs, including MK-1775 (adavosertib), paclitaxel, pevonedistat, ULK1_4989, vinblastine, and vinorelbine.Previous studies have demonstrated that circadian clockcontrol metabolism is a hallmark of cancer [12], and understanding the regulatory mechanisms of circadian rhythm links to the tumor immune microenvironment and metabolism could provide therapeutic benefts against tumors [63].

Conclusion
In conclusion, we integrated scRNA-seq data and bulk RNAseq data to comprehensively explore the immune characteristics in TME and identify the prognostic CRRGs used for survival prediction.Our data provided novel "timedependent" therapeutic options for HCC treatment.12 European Journal of Cancer Care Score Correlated with Poor Survival of HCC Patients.Te detailed fowchart of this study is shown in Figure S1.Among the TCGA-LIHC cohort, a total of 8831 DEGs (7793 upregulated and 1038 downregulated) were screened between HCCs and normal samples with the thresholds of |log 2 FC| > 0.585 and P value of <0.05 (Figures 1(a) and 1(b), and Table Figures 4(k)-4(m) and 4(p)-4(r)).Te Kaplan-Meier survival curves indicated that the high-risk group showed worse survival than the low-risk group (Figures 4(i), 4(n), and 4(s)).AUC values of the ROC curves for the 1-, 3-, and 5year survival analyses indicated the accuracy of the model for the survival prediction (Figures 4(j), 4(o), and 4(t)).

Figure 1 :
Figure 1: High-CR score correlated with poor survival of HCC patients.(a) Te volcano plot of DEGs between tumor and normal samples of HCCs in the TCGA-LIHC cohort.(b) Heatmap of DEGs between tumor and normal samples of HCCs in the TCGA-LIHC cohort.(c) Venn plot of the diferentially expressed CRRGs by overlapping DEGs and CRRGs.(d) Heatmap of 40 intersected CRRGs between tumor and normal samples of HCCs in the TCGA-LIHC cohort.(e) Kaplan-Meier plots for HCCs in high-CR score and low-CR score groups.

) and 6
(b)).Terefore, a nomogram was constructed by combining risk score and clinical characteristics (Figure 6(c)).Te calibration curves indicated the efciency of the nomogram for survival prediction (Figure 6(d)).Time-dependent ROC curves indicated the accuracy of the nomogram for survival prediction (Figure 6(e)).Te DCA curves revealed the discriminative ability of the nomogram for survival prediction (Figure 6(f )).

Figure 2 :
Figure 2: Identifcation of the key modules and hub genes related to circadian rhythm.(a) Sample clustering to identify the outliers.(b) Te scale-free index for soft-thresholding powers.(c) A dendrogram based on hierarchical clustering with optimal soft thresholds.(d) Heatmap of the correlation between modules and features (CR score).(e) Te histogram of the gene signifcance (GS) across modules.(f, g) Scatter plot of the correlation between module membership (MM) and gene signifcance (GS) of score in brown and green modules, respectively.

Figure 3 :
Figure 3: Dimensionality reduction, cell annotation, and identifcation of the CR score-related cell subpopulations.(a) Left: the UMAP of the 28 cell clusters.Right: the UMAP of 11 cell types was annotated by diferent markers.(b) Te bubble plot shows the GO-BP analysis of DEGs for each cell type.(c) Left: the UMAP of the CR score for each cell type.Right: the UMAP of the cell types distributed into the high-CR score and low-CR score groups.(d, e) Te bubble plot shows the GO (BP, CC, and MF) and KEGG enrichment analyses based on upregulated and downregulated DEGs between high-CR score and low-CR score groups, respectively.

Figure 4 :Figure 5 :
Figure 4: Construction and validation of a circadian rhythm-related signature.(a) Venn plot of the diferentially expressed CRRGs by overlapping DEGs from scRNA-seq and bulk RNA-seq data.(b) Te forest plot shows the prognostic CRRGs by univariate Cox analysis.(c) Te trajectory of each independent variable with lambda.(d) Plots of the coefcient distributions for the logarithmic (lambda) series for parameter selection (lambda).(e) Te selected prognostic genes and corresponding coefcients.(f, k, and p) Te survival status of patients and risk score of the TCGA-LIHC cohort, LIRI-JP cohort, and GSE14520 cohort.(g, l, and q) Te survival status of the patients with highrisk and low-risk score groups in the TCGA-LIHC cohort, LIRI-JP cohort, and GSE14520 cohort.(h, m, and r) Te heatmap of the fve CRRGs between high-risk and low-risk score groups in the TCGA-LIHC cohort, LIRI-JP cohort, and GSE14520 cohort.(i, n, and s) Kaplan-Meier plot for HCCs in the high-CR score and low-CR score groups of the TCGA-LIHC cohort, LIRI-JP cohort, and GSE14520 cohort.(j, o, and t) ROC curves for 1-, 3-, and 5-year survival prediction in the TCGA-LIHC cohort, LIRI-JP cohort, and GSE14520 cohort.

Figure 6 :
Figure 6: Development of a predictive nomogram.(a, b) Univariate and multivariate Cox analyses were performed to identify independent risk factors.(c) A nomogram for survival prediction based on risk factors and clinical characteristics.(d-f) Te calibration curve, DCA curve, and the time-dependent ROC curve were used to evaluate the efciency, predictive ability, and performance of the predictive model.

Figure 7 :Figure 8 :
Figure 7: Biological function analyses.(a-c) Bubble plots of the GO enrichment analysis, including BP, CC, and MF.(d) Bubble plots of the KEGG pathway enrichment analysis.(e) Bubble plots of the hallmark pathway enrichment analysis.

Figure 9 :
Figure 9: Prediction of the response to immunotherapy and chemotherapy in high-risk and low-risk groups.(a, b) Heatmap and histogram of the diferentially expressed immune checkpoint-related genes between high-risk and low-risk groups.(c, d) Histogram of the diferentially expressed TIDE value and IPS between high-risk and low-risk groups.(e, f ) Te diferential IC 50 values of chemotherapeutic drugs between high-risk and low-risk groups.