Integrated machine learning constructed a circadian-rhythm-related model to assess clinical outcomes and therapeutic advantages in hepatocellular carcinoma

Introduction

Hepatocellular carcinoma (HCC) is one of the most lethal malignancies around the world (1). Viral hepatitis, excessive alcohol intake, moldy food consumption, and metabolic liver disease all contribute to HCC (2). With increasingly expanded indications, surgical resection remains the mainstay of HCC treatment. Unfortunately, even after surgical intervention and active postoperative management, HCC still carries a dismal prognosis with 5-year recurrence rate of 70% (3). Therefore, identifying novel and reliable HCC biomarkers to improve patient outcomes, through accurate prognostic assessment and individualized treatment, is a critical but presently unmet clinical need.

Circadian rhythm (CR) is a ubiquitous time-keeping system coordinating physiological processes (wake-sleep, eating-fasting and activity-rest) with the daily cycles of daylight and darkness throughout a 24-hour period (4,5). At the molecular level, CR regulates cell cycle, DNA repair, apoptosis, proliferation, autophagy, metabolic nodes and immune components (6). Long-term irregular CR could increase the risks of a series of chronic diseases, especially cancers (7). According to epidemiological and experimental evidence, CR disruption has emerged as the hallmark of tumor onset and progression recently. For instance, chronic jet lag pushed fatty livers in wild-type mice, which ultimately developed to fibrosis and HCC (8). Moreover, the production of spontaneous circulating tumor cells (seeds of distant metastasis) coincides with sleep (9). These have raised hope for optimizing treatment (6). Wang et al. demonstrated that CR in dendritic cells (DCs) directly governed the anti-tumor response of CD8⁺ T cells in both mice and humans (10). However, the relationship between CR and cancer progression remains enigmatic, which may hamper recognizing underlying biomarkers to guide tailored therapies.

The last decade has witnessed a huge revolution in machine learning available for the entire spectrum of HCC management, encompassing risk assessment, diagnosis, and prognostic estimation (11). Considering the crucial role of machine learning in medical big data, we aimed to integrate machine learning into our research to excavate robust biomarkers to explore untapped potential correlation between CR and HCC progression. Here, HCC patients were divided into two CR-related phenotypes with distinct survival outcomes and functional annotations. The tumor microenvironment (TME) and treatment resistance of patients between two phenotypes were identified. Then, we screened the crucial circadian rhythm-related genes (CRRGs) and explored their prognostic values. Furthermore, the CR predictor was constructed employing machine learning algorithms to forecast the prognosis and therapeutic response. A nomogram was established based on the novel CR predictor; further, it may facilitate individualized survival prediction and pave the way for therapeutic regimens selection. Finally, we investigated the potential causal relationship between the model genes in CR predictor and HCC, using Mendelian randomization (MR). We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1155/rc).

Methods

Data acquisition

The RNA sequencing (RNA-seq) data and full clinical annotations of HCC patients were respectively extracted from The Cancer Genome Atlas (TCGA) Program and University of California Santa Cruz (UCSC) Xena platform. We also obtained expression matrices and corresponding clinical data from International Cancer Genome Consortium (ICGC) and the Gene Expression Omnibus (GEO) database. Gene list for CR was retrieved from the AmiGO 2 Web portal (http://amigo.geneontology.org/amigo) most recently. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Epigenetic and transcriptional features analysis

Based on the TCGA database, the types, rates, exclusive and co-occurrence of DNA mutations were explored via “maftools” package. Somatic copy number alteration (SCNA) data were downloaded from the cbioportal (https://www.cbioportal.org). Deep deletion and amplification were defined by values of −2 and 2 separately. The SCNA results and expression distributions were visualized by “ggplot2” and “ggpubr” package.

Unsupervised consensus clustering for CR regulators

Unsupervised Consensus clustering algorithm was employed to identify CR-mediated phenotypes via the “ConsensusClusterPlus” package. This clustering process was repeated for 1,000 times by resampling 80% items, categorizing each sample by k-means algorithm based on the Euclidean distance. To verify whether each phenotype was relatively independent, principal component analysis (PCA) was applied via “FactoMineR” R package. Using R package “survival” and “survminer”, the overall survival (OS) was appraised. Differentially expressed genes (DEGs) between two CR phenotypes was distinguished by “Limma” package. Furthermore, the Sankey diagram and heatmap were drawn to show the relationship between the DEGs expression, clinical characteristics, and CR phenotypes.

Functional enrichment and immune profile analysis

To probe into the variation of pathway activity in two CR phenotypes, Gene Set Variation Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA) were performed utilizing the “c2.cp.kegg.v2022.1.Hs.symbols.gmt” extracted from MSigDB database (12) (R package “GSVA” and “clusterProfiler”). Immune cell and stromal cell infiltration in HCC microenvironment were calculated by single-sample GSEA (ssGSEA) (13) and Microenvironment Cell Populations (MCP) counter (14) respectively. Beyond this, immune-related gene expression, mutant-allele tumor heterogeneity (MATH) values (15) and Tumor Immune Dysfunction and Exclusion (TIDE) score (16) were compared between two CR phenotypes. The 50% of the maximum inhibitory concentration (IC₅₀) for chemotherapeutic drugs were yielded by “oncoPredict” package.

Selection and verification of crucial CRRGs

The protein-protein interaction (PPI) network of CRRGs was established via Search Tool for the Retrieval of Interaction Gene (STRING) database with a reliable threshold (score >0.7). To elucidate key interacted genes, 4 algorithms of cytoHubba plug-in based on Cytoscape software were applied, including maximal clique centrality (MCC), maximum neighborhood component (MNC), edge percolated component (EPC) and Degree. Top 20 nodes (hub genes) scored by the 4 algorithms were overlapped by “UpSetR” R package. DEGs [identified using “Limma” R package with adjusted P value <0.001 and |log fold change (FC)| >1] in GSE208668 was used to validate hub genes. Ultimately, we intersected hub genes and DEGs to obtain crucial CRRGs. Relationship between expression and methylation of crucial CRRGs at pan-cancer level were assessed by Gene Set Cancer Analysis (GSCA) online tools (17). CpG methylation patterns was downloaded from MethSurv database (18). The integrated repository portal for tumor-immune system interaction database (TISIDB) enables us to look into expression relationship between CRY2 and chemokines or immune checkpoints through Spearman correlation coefficient (19). The distribution and expression of CRY2 in single-cell transcriptome (GSE140228 and GSE179795) were visualized by the Tumor Immunization Single Cell Hub (TISCH) database (20).

Cell culture and plasmid transfection

HCC cell lines (Huh7 and Hep3B) were obtained from Wuhan Pu-nuo-sai Life Technology Co. Ltd. (Wuhan, China). All the cell lines were cultured in Dulbecco’s modified Eagle’s medium containing 10% fetal bovine serum in a humidified atmosphere with 5% CO₂, at 37 ℃. The CRY2 overexpression plasmid (92539-1) and corresponding negative control vector plasmid (CON468) were synthesized by GeneChem (Shanghai, China). Lipofectamine 3000 (Invitrogen, USA) was applied to complete transfection.

Clone formation assay and wound healing assay

For the colon formation assay, cells were seeded in 6-well plates at a density of 1,000 per well. After 2 weeks, the cells were fixed in 4% paraformaldehyde and stained with 0.1% crystal violet. For the wound healing assay, cells were seeded in 6-well plates. At 90% confluence, the cell layer was scratched with a 200 µL tip. The migration was monitored and photographed at 0 and 24 h.

Cell cycle assay and apoptosis assay

For the cell cycle assay, cells were fixed with 70% pre-chilled ethanol overnight. Then, PI and RNase A working solution (Solarbio, China) were added. The rate of different cell cycle was detected by flow cytometry. For the apoptosis assay, cells were treated with fluorescein isothiocyanate (FITC)-Annexin V antibody and propidium iodide (PI) staining (Beyotime, China) in a dark environment. The data were analyzed via FlowJo (version 10.8.1).

Development and verification of the CR predictor via machine learning algorithms

HCC patients with complete survival information were enrolled, including TCGA cohort (n=363, training group), ICGC cohort (n=240, external validation group) and GSE14520 cohort (n=225, external validation group). Moreover, the expression of each cohort was normalized by a “scale” function. First, univariate Cox regression screened 66 CRRGs out for further analysis. Second, the random forest (RF) (21) and extreme gradient boosting (XGBoost) (22) algorithm were adopted to prioritize variables via “randomForest” and “xgboost” R packages. Briefly, the survival state for HCC patients was recognized as the response variable, and the CRRGs were selected as the explanatory variables. Majority voting was used to combine decision trees in RF algorithm. When an explanatory variable is added to the RF model, its importance is determined by the degree to which Gini coefficient falls. XGBoost algorithm is often regarded as an uninterpretable “black box”: the function between the features and the response is invisible to the researcher. Therefore, we leveraged SHapley Additive exPlanations (SHAP) (23), a game theoretic approach to explain how XGBoost worked. Thereafter, common candidate genes were selected from RF and XGBoost. Third, least absolute shrinkage and selection operator (LASSO) Cox regression shrunk candidates and minimized the risk of overfitting (R package “glmnet”, tenfold cross-validation and 1,000 re-sampling). Ultimately, CR predictor was calculated as follows:

Riskscore=∑i=19βi∗Ei[1]

where βi and Ei denoted the coefficient and expression of each gene, respectively.

Immunotherapeutic response and other drug sensitivity evaluation for CR predictor

First, Pearson correlation analysis screened the link between immune checkpoints and CR predictor. Then, a list of genes representing the seven-step anticancer immune process was obtained from the Tracking Tumor Immunophenotype (TIP) database (24), and enriched by GSVA algorithm. The R package “corrgram” linked risk scores with the genetic traits. The immune score and TME score were calculated by “ESTIMATE” and “TMEscore” package, respectively. The single-cell tumor immune microenvironment (scTIME) Portal, which houses single-cell RNA sequences (scRNA-seq) focusing on tumor immune microenvironment, was utilized to analyze model genes in GSE125449, GSE98638, and GSE140228. Further, we calculated risk score in GSE104580 and GSE109211 to probe into the capacity of CR predictor to forecast transcatheter arterial chemoembolization (TACE) or sorafenib efficacy. Last but not least, the Connectivity Map (CMap) web-based database was used to predict potential therapeutic agents and their underlying mechanisms (25).

Construction of a CR predictor-based nomogram

Nomogram is a pictorial representation of a scoring system for prognosis prediction. To excavate the independent prognostic potential of CR predictor in HCC, univariate and multivariate Cox analyses were conducted (R package “forestplot”). Age, stage and CR predictor integrated nomogram was generated by “cph” in “rms” package. The consistency index (C-index), calibration curve and decision curve analysis (DCA) were used to verify the predictive capacity for CR predictor (R package “rms” and “ggDCA”). Furthermore, to facilitate the use of nomogram, user-friendly dynamic nomogram was developed (R package “DynNom”).

MR analyses

To unveil the association between genes in CR predictor and the risk of HCC, MR analyses were performed. We obtained the Genome-Wide Association Studies (GWAS) data for HCC (finn-b-C3_LIVER_INTRAHEPATIC_BILE_DUCTS) from IEU Open GWAS and cis-expression quantitative trait locus (eQTL) data (26) for model genes. Instrumental variables closely associated with model genes (P<5×10⁻⁶) were selected as single nucleotide polymorphisms (SNPs) (r²<0.01 and distance >5,000 kb). Considering the less stringent threshold, F statistics were calculated. “TwoSampleMR” and “Coloc” R package was employed for causality and colocalization, respectively.

Statistical analysis

Data analyses were performed in R, GraphPad Prism software and FlowJo. The survival distribution was visualized by Kaplan-Meier (KM) plotter via the log-rank test. When comparing data from two groups, Wilcoxon test was applied for continuously distributed variables (R function “wilcox.test”), while Chi-square test was utilized for categorical variables (R function “chisq.test”). Correlation between two continuous variables was measured by Pearson’s correlation coefficient. For all analyses, P value <0.05 was deemed valuable.

Results

Epigenetic features of CR regulators in liver hepatocellular carcinoma (LIHC)

To explore the epigenetic features of CR regulators in HCC, the prevalence of somatic mutations of 212 CRRGs in LIHC was shown on Figure 1A. CR regulators mutated in 191 of the 258 (74.03%) LIHC samples, predominantly by missense mutations and frame shift deletion mutations. Among them, TP53 had the highest mutation rate, while PRKDC (7%), KMT2A (5%) and MYCBP2 (5%) followed (Figure S1A-S1D). Additionally, we investigated the mutation co-occurrence across CR regulators. The co-occurrence of ADCY1 and PRKDC, along with HNF4A and ZFHX3, were revealed (Figure 1B). Figure 1C displayed the copy number alteration (CNA) frequency of CR regulators in LIHC. CHRNB2, RORC, CIART, NTRK1 and FBXL6 exhibited widespread copy number amplification, whereas almost no deep deletion was exhibited in CRRGs with high CNA frequency. Moreover, compared to 50 normal tissues, CHRNB2, CIART, FBXL6, NCOA2, and HNRNPU expressed considerably higher in 369 liver cancer tissues (P<0.001, Figure 1D).

Figure 1 Genetic alterations and transcriptional expression of CRRGs in HCC. (A) The top 20 most frequently mutated CRRGs in TCGA. (B) Co-occurrence and mutually exclusive analyses of mutant CRRGs. (C) CNA frequency of circadian rhythm regulators in TCGA. blue dots, deep deletions; red dots, amplifications. (D) Expression distributions of high CNA CRRGs between tumor and normal samples in TCGA. *, P<0.05; **, P<0.01; ***, P<0.001; ns, no statistical significance (Wilcoxon test). CRRGs, circadian rhythm-related genes; HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; CNA, copy number alteration; TMB, tumor mutation burden.

Recognition of CR-related phenotypes in LIHC

Among the clustering variable (k) values, k=2 generated the best outcomes according to the clustering heatmaps and cumulative distribution function (CDF) curves (Figure 2A, Figure S2A,S2B). HCC patients were well categorized into Cluster 1 (141 cases, 38.21%) and Cluster 2 (228 cases, 61.79%). Similar to the thermograms, the distinct clusters were corroborated by PCA (Figure 2B). Besides, Figure 2C indicated that Cluster 2’s survival probability was higher than Cluster 1. The GSE14520 cohort also demonstrated significant differential survival outcomes (Figure S2C). The clinicopathological parameters and CRRG expression were also compared. As depicted in Figure 2D, Cluster 1 exhibited a larger prevalence of advanced histologic grade (G3 or G4), which was closely tied to poor OS status, as well as Child-Pugh classification grade and tumor node metastasis (TNM) stage. Most CR regulators, such as SFPQ, HNRNPL, HNRNPU, HDAC2, TOP2A, EZH2, CDK1 and PPP1CC were noticeably high-expressed in Cluster 1, while ASS1, KLF9, MTTP, F7 and RORC were the opposite (Figure 2E). The observed consistency implied the accuracy of the CR patterns in predicting the prognosis of HCC.

Figure 2 Phenotypes of circadian rhythm and their prognostic value in HCC divided by consistent clustering. (A) Two phenotypes of circadian rhythm were identified by k-means using 1,000 iterations in TCGA. (B) PCA to verify the two circadian rhythm-related phenotypes in TCGA. (C) The KM curve of circadian rhythm-related phenotypes in TCGA via the log-rank test. (D) The Sankey diagram portrayed the correlations across the circadian rhythm-related phenotype clusters, histologic grade, Child-Pugh classification grade, stage and survival status in TCGA. (E) CRRGs expression in two clusters and associated clinical characteristics. Differentially expressed CRRGs were identified by “limma” package with the threshold of adjusted P value <0.001. G1, Grade 1 (histologic grade); G2, Grade 2 (histologic grade); G3, Grade 3 (histologic grade); G4, Grade 4 (histologic grade); HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; CRRGs, circadian rhythm-related genes; PCA, principal component analysis; KM, Kaplan-Meier.

Functional comparison disclosing the link between CR phenotypes and immune landscape

To elucidate the discrepancies between enrichment pathways and CR clusters, GSVA enrichment was conducted (Figure 3A). Cluster 1 was enriched in immune and carcinogenic pathways, such as MAPK, ErbB, p53, mTOR, and Wnt signaling pathway. Cluster 2, on the other hand, was enriched in metabolism-related pathways. Consistent with the results from GSVA, GSEA analysis suggested that immune-related pathways were significantly activated in Cluster 1 (Figure 3B). So, we further evaluated the infiltration level of 28 immune cell types. Intriguingly, Cluster 1 showed abundant infiltration not only in immune-promoting cells, such as CD4 T cell, central memory CD8 T cell, immature B cell, memory B cell, Tfh, Th17 and Th2, but also in immune-inhibiting cells, like myeloid-derived suppressor cells (MDSC) and regulatory T cells (Treg). CD 56 bright natural killer (NK) cell and eosinophil, in contrast, were more prevalent in Cluster 2 (Figure 3C). Meanwhile, preferential expression of human leukocyte antigen (HLA) and major histocompatibility complex (MHC) gene sets were found in Cluster 1 (Figure S2D,S2E).

Figure 3 Functional enrichment and immune landscape of circadian rhythm-related phenotypes. (A) GSVA enrichment analysis of KEGG biological pathways for circadian rhythm-related phenotypes. Orange: activated pathways; blue: inhibited pathways. (B) The KEGG signaling pathways enriched by GSEA. (C) Immune cell infiltration between two phenotypes. *, P<0.05; **, P<0.01; ***, P<0.001; ns, no statistical significance (Wilcoxon test). GSVA, gene set variation analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, gene set enrichment analysis; NES, normalized enrichment score; PD-L1, programmed cell death-ligand 1; MDSC, myeloid-derived suppressor cells; ssGSEA, single-sample gene set enrichment analysis.

According to prior research, higher infiltration of immune cells was linked to the prognosis advantage (27). However, Cluster 1 with the higher immune cell infiltration did not exhibit a superior prognosis over Cluster 2. Current study suggested that interaction of cancer-associated fibroblasts (CAFs) with immune components could upregulate the expression of immune checkpoints, inducing T-cells dysfunction and evading immune surveillance (28). Therefore, we hypothesized that the antitumor immunity in Cluster 1 could be influenced by such factors. MATH score, a straightforward and quantitative reflection of intratumor heterogeneity, was found to be notably higher in the Cluster 1 (Figure 4A). Subsequent analysis confirmed that Cluster 1 was certainly featured with a larger fibroblast infiltration (Figure 4B). The variation in immune checkpoint expression across two clusters was also examined. Cluster 1 showed higher expression levels of PD-1, PD-L1, CTLA-4, CD86, TIGIT, CD47, TIM-3 and LAG-3 (Figure 4C), indicating a tendency for immune evasion. Meanwhile, Cluster 1 exhibited a higher TIDE score than Cluster 2, representing a greater possibility for immune escape and fewer benefits from immunotherapy (Figure 4D,4E). Moreover, we utilized GDSC2 database to assess the IC₅₀ of drugs, covering sorafenib, 5-fluorouracil, cisplatin, vinblastine, cyclophosphamide and irinotecan. IC₅₀ values of these drugs were lower in Cluster 1, implying that Cluster 1 were more sensitive to chemotherapeutic and targeted agents (Figure 4F).

Figure 4 Prediction of circadian rhythm-related phenotypes in treatment benefits. (A) MATH predicting intratumor heterogeneity between two phenotypes. (B) Infiltration of fibroblasts between two phenotypes. (C) Expression of immune checkpoints between two phenotypes (Wilcoxon test). (D,E) The distribution of TIDE scores predicting response of immunotherapy between two phenotypes. (F) Response to six common chemotherapeutic and targeted drugs between two phenotypes. The y-axis represents the IC₅₀ for each drug. The differences were compared by Wilcoxon test (violin chart) and Chi-squared test (percent stacked bar chart). **, P<0.01; ***, P<0.001 (Wilcoxon test). MATH, mutant-allele tumor heterogeneity; TIDE, Tumor Immune Dysfunction and Exclusion; IC₅₀, 50% of the maximum inhibitory concentration.

Screening of crucial CRRGs

A PPI network (211 nodes, 1,680 edges) was imported into Cytoscape for further analysis (Figure 5A). Based on the intersection of the top 20 genes through 4 topological algorithms (MCC, MNC, EPC and Degree), 11 hub genes were obtained (BTRC, PER2, NR1D1, NPAS2, CRY1, ARNTL, PER1, CRY2, CLOCK, CSNK1E, EP300) (Figure 5B). However, the results based solely on the algorithm might be weak in real-world biological relevance and trustworthiness. Desynchrony between the internal (organisms) and external (the earth) CR will bring on insomnia (29). Thus, we identified 1,348 DEGs between older people with and without insomnia disorder from GSE208668 (Figure 5C). Ultimately, by integrating the DEGs and 11 hub genes, three crucial CRRGs (CRY1, CRY2 and EP300) were detected (Figure 5D). The expression of crucial CRRGs was inversely linked with their methylation levels in multiple cancers (Figure 5E).

Figure 5 Recognition of crucial CRRGs. (A) The PPI network for all circadian rhythm regulators. (B) The identification of hub genes through intersecting across top 20 nodes scored 4 topological algorithms (MCC, MNC, EPC and Degree). The height of each bar represented the number of genes in that overlap. Eleven common hub genes from the first bar were extracted. (C) Volcano plots showing DEGs in GSE208668. (D) The overlaps of 11 hub genes and DEGs. (E) Relationship between expression and methylation of three crucial CRRGs in pan-cancer. CRRGs, circadian rhythm-related genes; PPI, protein-protein interaction; MCC, maximal clique centrality; MNC, maximum neighborhood component; EPC, edge percolated component; COL, color; DEGs, differentially expressed genes; FDR, false discovery rate; mRNA, messenger RNA.

To further verify the prognostic value of crucial CRRGs, univariate Cox regression was conducted in TCGA. Only CRY2 served as an effective protective variable [hazard ratio (HR) =0.714, 95% confidence interval (CI): 0.568–0.896, and P=0.004, Figure 6A]. According to the KM curve, higher expression of CRY2 in TCGA (logrank P<0.001, Figure 6B) and ICGC (logrank P=0.005, Figure 6C) cohorts was associated with superior OS, suggesting that CRY2 expression may impede HCC development. Moreover, to investigate the independence of CRY2, multivariate Cox regression was then performed (Figure 6D). After adjustment for pathologic T and histologic grade, CRY2 still exhibited robust prognostic assessment ability (HR: 0.741, 95% CI: 0.591–0.929, P=0.01). Then we performed the relevant basic assay validation for CRY2. Overexpression of CRY2 inhibited clone formation capability in Huh7 (P=0.02) and Hep3B cells lines (P<0.001, Figure 6E). In the wound-healing assay, we found that overexpression of CRY2 inhibited wound closure speed in both Huh7 (P<0.001, Figure 6F) and Hep3B cells (P<0.001, Figure 6G). Employing Annexin V and PI staining for flow cytometry, we observed that CRY2 overexpression induced apoptotic responses (Figure 6H). HCC cells with overexpressed CRY2 exhibited an inhibitory effect on cell cycle progression (Figure 6I). Therefore, we concentrated on the potential mechanism by which CRY2 predicts prognosis in HCC. First, CRY2 methylation map revealed 16 CpG sites, three of which displayed impressive correlations with prognosis, including cg02208899 (P=0.01), cg21002528 (P=0.045), and cg22234194 (P=0.03, Figure S3A,S3B). Next, we discovered that CRY2 was negatively correlated with immune inhibitors (CTLA4, HAVCR2, PDCD1 and TIGIT), and was positively associated with chemokines (CCL14, CCL15, CCL16 and CCL25) (Figure S4A,S4B). On top of that, we annotated the GSE140228 and GSE179795 single cell expression based on the TISCH2 database. In GSE140228 dataset, CRY2 was chiefly expressed in conventional CD4 T cells, proliferating T cells, exhausted CD8 T cells, NK, innate lymphoid cells (ILC), DC, B, plasma, mono/macro and mast cells (Figure S4C). In GSE179795 dataset, CRY2 was predominantly enriched in CD8 T cells, NK, ILC, and mono/macro (Figure S4D).

Figure 6 Prognostic value for crucial CRRGs. (A) Univariate Cox analysis of three crucial CRRGs with OS in TCGA. (B,C) Kaplan-Meier curves for CRY2 in TCGA and ICGC. The median expression of CRY2 was set as the cutoff point. (D) Multivariate Cox analysis of variables including CRY2 expression, pathologic T and histologic grade with OS in TCGA. The upper part represents hazard ratio, and the lower part in parentheses represents lower 95% confidence interval and upper 95% confidence interval. The x-axis represents the hazard ratios associated with HCC survival in multivariate Cox analysis. (E-G) Overexpression of CRY2 inhibited the proliferation (E) and migration abilities (F,G) of Huh7 and Hep3B cells. (E) is stained by crystal violet. The magnification of (E) is 1×. The magnification of (F,G) is 100×. The observation method of (F,G) is microscope. (H,I) Overexpression of CRY2 induced apoptotic responses (H) and inhibit cell cycle (I) in Huh7 and Hep3B cells. *, P<0.05; **, P<0.01; ***, P<0.001. CRRGs, circadian rhythm-related genes; OS, overall survival; HR, hazard ratio; CI, confidence interval; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; PI, propidium iodide; G1, first gap; G2, second gap; S, synthesis; RCS, restricted cubic spline.

Construction and validation of CR predictor

To establish a clinically applicable model of the CR related signature for prognosis, we employed three machine learning algorithms on 363 TCGA-LIHC patients with survival information. The study flowchart was illustrated in Figure 7A. After univariate Cox regression analysis, 66 prognosis-related CRRGs were selected for further analysis (Figure S5). Subsequently, RF algorithm was employed to identify top 20 crucial genes. The relative importance of variable was quantified by the mean decrease in the Gini index (Figure 7B). The higher mean decrease in Gini index indicates higher variable importance. Next, we calculated the SHAP value of XGBoost algorithm (Figure 7C). The higher a feature’s SHAP value, the greater its prognostic prediction potential. The SHAP dependency analysis offered an insightful method for interpreting how a single feature influences the output of the so-called ‘black box’ XGBoost prediction model (Figure S6). Finally, the top 20 genes selected by the RF and XGBoost algorithm separately were overlapped, yielding 14 genes as survival-related potential candidate biomarkers (Figure 7D). Ultimately, we narrowed down the 14 genes implementing LASSO Cox regression algorithm with tenfold cross-validation (Figure 7E,7F). Nine genes with non-zero coefficients were incorporated into the predict model with minimum lambda (0.029). The computational equation of CR predictor was: Risk Score = (0.06110786 * SFPQ) + (0.284424475 * PPP1CB) + (0.238875693 * NR1H3) + (−0.17090971 * CRY2) + (−0.104810559 * PPARGC1A) + (0.085020819 * HDAC2) + (0.039030674 * PRKAA2) + (0.044803425 * HOMER1) + (0.050396753 * TNFRSF11A) (Figure 7G). Unified threshold to divide patients was established using the median. As anticipated, patients in Cluster 1 with poor prognoses primarily belonged to high-risk group (Figure 7H).

Figure 7 Construction of the CR predictor. (A) The workflow of the development of CR predictor. (B) RF algorithm and the top 20 important genes ranked by mean decrease Gini. When a new gene (explanatory variable) is added to the RF model, its importance to predict the survival status (response variable) is measured by the degree to which the Gini coefficient falls. The length of the line segment and the size of the point represent the value of the mean decrease Gini. (C) XGBoost algorithm and the top 20 important genes interpreted by SHAP values. The value next to the y-axis was the mean SHAP value. The higher the SHAP value of a gene (explanatory variable) was given, the higher risk of death (response variable) the patient would have. Each patient was allocated one dot in each row. The x-position of the dot depended on its SHAP value (output prediction value); the color of the dot depended on each patient’s gene expression (input observation value). Purple and yellow dots corresponded to higher and lower gene expression level respectively. (D) Venn plot of common candidate biomarkers from RF and XGBoost algorithm. (E,F) LASSO Cox regression analysis of 14 common candidate biomarkers by tenfold cross-validation and 1,000 re-sampling. (G) Corresponding coefficients for the 9 genes that make up the CR predictor. (H) The relationships among the circadian rhythm-related phenotype clusters, CR predictor and survival status. CR, circadian rhythm; RF, random forest; XGBoost, extreme gradient boosting; SHAP, SHapley Additive exPlanations; LASSO, least absolute shrinkage and selection operator; Coef, coefficient.

To further validate CR predictor, independent HCC cohorts were enrolled. As the risk score increased, the probability of mortality in HCC patients rose. Besides, the expression of nine model genes was comparable (Figure 8A,8B). Patients with higher risk score suffered from worse prognosis significantly (P<0.001) in both TCGA datasets (Figure 8C) and two external validation cohorts (Figure 8D, Figure S7A). The TCGA cohort’s 1-, 3-, and 5-year survival area under the curve (AUC) reached 0.762, 0.746, and 0.746, respectively (Figure 8E), and they were 0.644, 0.619, and 0.637 for GSE14520 (Figure S7B). Analogously, 1-, 2-, and 3-year survival AUC for ICGC were 0.686, 0.698, and 0.686 (Figure 8F), suggesting the CR predictor was reliable for forecasting progression. CR predictor was notably related to HCC prognosis, even after adjusting for other variables (Figure S7C,S7D). These all disclosed the CR predictor could recognize HCC patients with favorable prognosis.

Figure 8 Survival evaluation of CR predictor in the training dataset (TCGA) and the external validation dataset (ICGC). (A,B) Distribution and expression of 9 model genes. (C,D) Kaplan-Meier curves for CR predictor (log-rank test). (E,F) Time-dependent ROC curves of CR predictor. CR, circadian rhythm; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; ROC, receiver operating characteristic; AUC, area under the curve.

CR predictor-based therapeutic regimens for HCC patients

Immunotherapy has revolutionized the management of HCC (30). Pearson’s correlation revealed a substantial positive connection between risk score and immune checkpoints (PD-1, PD-L1, CTLA-4, CD86, TIGIT, LAG-3, CD47, TIM-3) in TCGA and ICGC datasets (Figure 9A). To further elucidate the underlying immunotherapeutic potential, we calculated enrichment scores for seven-step cancer-immunity cycle (Figure 9B). Strikingly, risk score was positively related to Steps 1, 5 and 6. The high-risk group also exhibited higher immune and TME score, indicating the complexed tumor immune microenvironment (Figure 9C,9D). We then explored detailed distribution of model genes in HCC patients by scRNA-seq data (GSE125449, GSE98638 and GSE140228). Consistently, PPP1CB, HDAC2 and SFPQ, owning positive coefficients in CR predictor, expressed highly in various immune cells (Figure S8). Altogether, CR predictor closely correlated with the tumor immune microenvironment, and may count for predicting HCC immune response.

Figure 9 The immunotherapeutic potential for CR predictor. (A) Pearson correlation between CR predictor and immune checkpoints expression. The width of the inside ribbons denoted the Pearson correlation coefficients (R); the color of the outside tracks denoted the P value. (B) Correlation between CR predictor and seven-step cancer-immunity cycle in TCGA. The Pearson correlation coefficients (R) were rendered in the wedge-like sectors of the pie charts. The red color represented a positive correlation. The deeper the red, the bigger the R value. (C) Distribution of immune scores in TCGA (Wilcoxon test). (D) Distribution of TME scores in TCGA (Wilcoxon test). CR, circadian rhythm; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; TME, tumor microenvironment.

TACE and sorafenib are recommended for HCC patients (31,32). We investigated whether CR predictor could identify HCC patients’ response to chemotherapy and targeted treatment. The GSE104580 cohort of HCC patients undergoing TACE therapy was examined (Figure 10A,10B), and high-risk patients showed a higher tendency for TACE resistance (P<0.001). The AUC of the CR predictor in forecasting TACE is 0.688 (Figure 10C). In GSE109211, patients who response to sorafenib exhibited a lower risk score (Figure 10D). Low-risk patients were more sensitive to sorafenib (P<0.001, Figure 10E). The receiver operating characteristic (ROC) curves further confirmed the superior performance of the CR predictor in forecasting sorafenib response (AUC =0.877, Figure 10F).

Figure 10 Drug sensitivity and underlying compounds targeting the CR predictor. (A,D) Difference in risk score among (A) TACE and (D) Sorafenib responder and non-responder (Wilcoxon test). (B,E) The proportion of HCC patients with distinct responses to (B) TACE and (E) Sorafenib therapy (Chi-squared test). (C,F) The ROC curves of the CR predictor for forecasting (C) TACE and (F) Sorafenib response. The data in (C,F) indicate the best cut-off point value and its specificity and sensitivity for the ROC curve. (G) Potential small-molecule compounds and their drug mechanisms. X-axis showed drug names, and y-axis showed mechanisms targeted by drugs. Forty-one compounds regulating HCC cell lines were sorted from the CMap database. CR, circadian rhythm; TACE, transcatheter arterial chemoembolization; ROC, receiver operating characteristic; HCC, hepatocellular carcinoma; AUC, area under the curve; CMap, Connectivity Map.

The CMap database offers insight into intricate connections between gene expression, phenotype, compounds and mode of action (MoA). So, we intended to find novel compound candidates that could target the CR predictor via the CMap database. The top 150 upregulated and 150 downregulated DEGs were obtained from two risk groups. With the criteria of normalized connectivity score >1.4, we identified 41 compounds regulating HCC cell lines (Huh7 and HepG2). Figure 10G revealed that 5 compounds (albendazole, D-64131, mebendazole, nocodazole and oxibendazole) shared the mechanism of tubulin inhibitor. Most drug-target MOA belong to the inhibitory category. For example, acemetacin and tolmetin shared the mechanism of cyclooxygenase inhibitor; bendroflumethiazide and hexamethylene amiloride shared the mechanism of sodium channel inhibitor.

Development of CR predictor-based nomogram model for individualized evaluation

A nomogram incorporating risk score, age and stage was established with a C-index of 0.723 (Figure 11A). The nomogram also exhibited satisfactory agreement between observed and predicted values (Figure 11B). Meanwhile, DCA verified the net benefit of the nomogram in a clinical context (Figure 11C). To quantify the Nomogram model’s accuracy, ROC curves were drawn in three independent cohorts (Figure 11D-11F). Based on the evidence above, our model served as a reliable and accurate prognostic indicator in clinical decision-making. An online survival calculator was developed for users to explore (https://lihc.shinyapps.io/CR_predictor/). The calculator collected input data including risk score, age and stage, and then automatically exported estimated survival profile for an HCC patient.

Figure 11 Construction of the nomogram integrating CR predictor and clinical parameters. (A) Nomogram for predicting OS in HCC. (B) Calibration curves of the nomogram. (C) DCA to evaluate nomogram’s clinical utility. (D-F) ROC curves of the nomogram in TCGA, ICGC and GSE14520. CR, circadian rhythm; HCC, hepatocellular carcinoma; OS, overall survival; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; DCA, decision curve analysis; ROC, receiver operating characteristic; AUC, area under the curve.

Causal effects of genes in CR predictor on HCC via MR analyses

To probe into the potential causal effect between model genes and HCC, MR analyses were conducted. Six cis-eQTLs were associated with the expression of 9 model genes after clumping SNPs (threshold: P<5×10⁻⁶, distance >5,000 kb, r²<0.01). Mean F statistics of the SNPs ranged from 21.04 to 2,545.35, indicating that the instrumental variables were powerful (Table S1). Among 9 model genes, only PPARGC1A showed a significant causal effect on HCC (Figure 12A-12D, Table S2). In the sensitivity analysis, no heterogeneity (P=0.66) and pleiotropy (P=0.30) were found (Table S3). A case-control study has also confirmed that PPARGC1A increased the risk of HCC in Moroccan population (33). Reverse causality (Figure S9A) and colocalization analyses (PP.H4 =4.56e−02, Figure S9B) were not supported.

Figure 12 Causal associations of model genes in CR predictor on HCC via MR analyses. (A) The exposure was cis-eQTLs for 9 model genes, and the outcome was GWAS for HCC; 6 cis-eQTLs were associated with the expression of 9 model genes after clumping SNPs. IVW method was selected. Only PPARGC1A showed a significant causal effect on HCC. (B) Scatter plot showing the association of the SNP effects on model genes in CR predictor against the SNP effects on HCC. (C) Funnel plot to visually assess heterogeneity. (D) Causal effects of PPARGC1A on HCC via reverse MR analyses. CR, circadian rhythm; HCC, hepatocellular carcinoma; MR, Mendelian randomization; SE, standard error; GWAS, Genome-Wide Association Studies; eQTL, expression quantitative trait locus; SNP, single nucleotide polymorphism; nSNP, number of single nucleotide polymorphism; OR, odds ratio; CI, confidence interval; IVW, inverse-variance weighted.

Discussion

As one of the most aggressive malignancies, HCC is contributing to the public health burden exponentially (34). Currently, several adjuvant therapeutic strategies have been approved for HCC, but their benefits are modest (35). Therapeutic breakthroughs are still warranted. CR is a global regulatory system that controls several cell processes, such as DNA damage response, redox homeostasis, and replicative immortality (36,37). The central circadian clock synchronizes liver clocks and controls the cyclic expression of transcription factors, different enzymes and metabolic master regulators (38,39). Therapeutic and lifestyle interventions that support a functional circadian clock can have significant health benefits. Limiting daily eating to a 6- to 10-hour window aligned with the body’s internal clock can enhance metabolic health and combat metabolic disorders (40). Exercising at different times of the day could be a beneficial strategy for managing metabolic diseases, as it helps align the body’s clock and regulate metabolism (41). When it comes to the therapeutic interventions, circadian clock-regulating biomaterials can promote disc regeneration by regulating the circadian clock of nucleus pulposus cells (42). Moreover, efficacy of immunotherapy can be improved by scheduling the treatment at specific times of the day (43). However, although molecular crosstalk between CR and HCC was emphasized, the deep-seated relationship between CRRGs and HCC patients’ prognosis has not been well illustrated.

In the present research, we extracted the CRRGs of HCC patients through unsupervised representation learning strategy, and defined two CR phenotypes (Cluster 1 and Cluster 2) with distinguished survival outcomes. The two phenotypes also showed heterogeneity in clinical characteristics, functional enrichment, TME characteristics, and drug sensitivity. Cluster 2 with longer survival showed an enrichment of metabolic pathways. Metabolic control by the CR underpinned cancer metabolism (44). Consistent with our results, feeding-driven fatty acid metabolism in liver, which was regulated by hepatic acetyl-CoA carboxylase (ACC), followed a 24-hour rhythm (45). Cluster 1 with poor prognosis was characterized by tumor-promoting signals. The activation of mTOR signaling pathway, which was also rhythmic (sensed by food intake), played a predominant role in hepatocarcinogenesis (46). These may explain the different tumor progression of HCC patients. Interestingly, Cluster 1 not only had a tight connection with effector immune cell receptor pathways (immune activation), but also presented a higher fibroblasts infiltration and immune checkpoints expression (immune suppression). The paradoxical results implied the high cellular heterogeneity in HCC (47). We speculated that the Cluster 1 patients retain their antitumor immunity potentiality, which might be blocked by the immunosuppressive ecosystem (48,49). According to our research, Cluster 1 patients may derive greater benefit from chemotherapeutic and targeted drugs. Previous reports have suggested that chronotherapy (an approach based on the CR) could minimize cisplatin-induced renal injury (50). CR modulation of 5-fluorouracil’s infusion rate may further enhance its therapeutic efficacy (51). Together, we confirmed the existence of CR-related phenotypes in HCC and offered a promising direction for clinical management.

Subsequently, with stringent selection criteria, we finally obtained three essential CR-related genes (CRY1, CRY2, and EP300). Cryptochrome 1 (CRY1), a nuclear coregulator, promoted tumorigenesis by regulating DNA repair and the G2/M transition temporally (52). HCC cells were profoundly vulnerable to perturbation of E1A binding protein P300 (EP300), which could elevate AFP protein level in plasma and promote metastasis in HCC patients (53). Cryptochrome 2 (CRY2) was the only protective factor in HCC among three crucial genes. In concordance with our results, CRY2 displayed anti-tumorigenesis effects in several cancers (54,55). Missense mutations of mouse CRY2 robustly repressed the growth-promoting effect of p53 (56). Moreover, ablation of CRY2 methylation inhibited cancer cell proliferation and invasion in vitro (57). These studies suggested that CRY2 served as a pivotal tumor suppressor.

Considering the impact of CR-related phenotypes and CRRGs on the clinical outcomes, three machine learning algorithms (RF, XGBoost and LASSO) were employed, and CR predictor was established. CR predictor comprised nine genes, including PPP1CB, NR1H3, HDAC2, SFPQ, TNFRSF11A, HOMER1, PRKAA2, PPARGC1A and CRY2, which were related to cancer initiation, and progression metastasis. PPP1CB controlled nuclear shape of cancer cells; the deletion of PPP1CB broke the nuclear envelope and triggered genome instability (58). NR1H3 promoted HCC progression via an HNF4α-dependent reciprocal network (59). HDAC2 facilitated pancreatic cancer metastasis by regulating TGFβ pathway (60). SFPQ played an important role in controlling the oncogenic splicing switch in HCC (61). TNFRSF11A regulated Ca²⁺-calcineurin/NFATC1-ACP5 axis to accelerate colorectal cancer cell metastasis (62). Expression of HOMER1 downregulated in HCC and was linked to tumor size (63). PRKAA2 was closely associated with 5-fluorouracil chemoresistance in gastric cancer (64). PPARGC1A induced mitochondrial biogenesis and oxidative phosphorylation to promote distant metastases in cancer cells (65). We also evaluated whether CR predictor was connected with immunomodulators, microenvironment, and treatment response. CR predictor can not only monitor the clinical outcomes for HCC, but also guide treatment. Finally, the nomogram, incorporating CR predictor, age and stage, exhibited accurately predictive performance for OS of HCC patients. Our nomogram is easy to use, and it may have broad development prospects for individualized survival prediction and therapeutic regimens selection in HCC patients.

However, limitations still exist. First, the prognosis of HCC patients varies depending on the etiology (viral hepatitis, excessive alcohol intake, moldy food consumption, and metabolic liver disease). Due to data restriction, we cannot take account of etiology’s impact on the HCC prognosis. Second, considering that the sequencing data were obtained from liver tissue, CR predictor may not apply to peripheral blood samples. Third, the crucial CRRGs identified in our study have not been verified in molecular experiments, which are meaningful works in the future. Last but not least, while our signature has been rigorously tested across numerous training and validation datasets to assess its predictive capabilities, it’s important to note that all the cohorts analyzed were retrospective. Our model should be further confirmed in prospective studies with large sample sizes.

Conclusions

To sum up, we depicted the stratification value of CR regulators in HCC and screened out the pivotal CRRGs. Employing machine learning algorithms, we developed a survival risk classifier named CR predictor, which could assess clinical outcomes and drug sensitivity for HCC patients accurately. Targeting CR could be a promising strategy against HCC.

Acknowledgments

The authors gratefully acknowledge contributions from the TCGA, ICGC and GEO databases.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1155/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1155/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1155/prf

Funding: This work was supported by the State Administration of Chinese Medicine Project (No. 20085-9-1), Project initiated by researchers from Jiangsu Province Hospital of Chinese Medicine (No. ZDXYS202208-1), Jiangsu Province Hospital of Chinese Medicine Peak Academic Talent Project (No. y2021rc19), Jiangsu Provincial Health and Medical Committee Key Projects (No. ZD2022070), Science and Technology Project of Affiliated Hospital of Nanjing University of Chinese Medicine (No. Y2020CX62), State Administration of Chinese Medicine Project (No. 20085-9-3), and Jiangsu Provincial Science and Technology Department Project (No. BE2019771).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1155/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Yang JD, Hainaut P, Gores GJ, et al. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat Rev Gastroenterol Hepatol 2019;16:589-604. [Crossref] [PubMed]
3. Tabrizian P, Jibara G, Shrager B, et al. Recurrence of hepatocellular cancer after resection: patterns, treatments, and prognosis. Ann Surg 2015;261:947-55. [Crossref] [PubMed]
4. Zhou L, Zhang Z, Nice E, et al. Circadian rhythms and cancers: the intrinsic links and therapeutic potentials. J Hematol Oncol 2022;15:21. [Crossref] [PubMed]
5. Halberg F. Physiologic 24-hour periodicity; general and procedural considerations with reference to the adrenal cycle. Int Z Vitaminforsch Beih 1959;10:225-96. [PubMed]
6. Sulli G, Lam MTY, Panda S. Interplay between Circadian Clock and Cancer: New Frontiers for Cancer Treatment. Trends Cancer 2019;5:475-94. [Crossref] [PubMed]
7. El-Athman R, Relógio A. Escaping Circadian Regulation: An Emerging Hallmark of Cancer? Cell Syst 2018;6:266-7. [Crossref] [PubMed]
8. Kettner NM, Voicu H, Finegold MJ, et al. Circadian Homeostasis of Liver Metabolism Suppresses Hepatocarcinogenesis. Cancer Cell 2016;30:909-24. [Crossref] [PubMed]
9. Diamantopoulou Z, Castro-Giner F, Schwab FD, et al. The metastatic spread of breast cancer accelerates during sleep. Nature 2022;607:156-62. [Crossref] [PubMed]
10. Wang C, Barnoud C, Kizil B, et al. Dendritic cells direct circadian anti-tumor immune responses. Ann Oncol 2022;33:S560. [Crossref]
11. Calderaro J, Seraphin TP, Luedde T, et al. Artificial intelligence for the prevention and clinical management of hepatocellular carcinoma. J Hepatol 2022;76:1348-61. [Crossref] [PubMed]
12. Liberzon A, Birger C, Thorvaldsdóttir H, et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015;1:417-25. [Crossref] [PubMed]
13. Barbie DA, Tamayo P, Boehm JS, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 2009;462:108-12. [Crossref] [PubMed]
14. Becht E, Giraldo NA, Lacroix L, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 2016;17:218. Erratum in: Genome Biol 2016;17:249. [Crossref] [PubMed]
15. Mroz EA, Rocco JW. MATH, a novel measure of intratumor genetic heterogeneity, is high in poor-outcome classes of head and neck squamous cell carcinoma. Oral Oncol 2013;49:211-5. [Crossref] [PubMed]
16. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
17. Liu CJ, Hu FF, Xie GY, et al. GSCA: an integrated platform for gene set cancer analysis at genomic, pharmacogenomic and immunogenomic levels. Brief Bioinform 2023;24:bbac558. [Crossref] [PubMed]
18. Modhukur V, Iljasenko T, Metsalu T, et al. MethSurv: a web tool to perform multivariable survival analysis using DNA methylation data. Epigenomics 2018;10:277-88. [Crossref] [PubMed]
19. Ru B, Wong CN, Tong Y, et al. TISIDB: an integrated repository portal for tumor-immune system interactions. Bioinformatics 2019;35:4200-2. [Crossref] [PubMed]
20. Sun D, Wang J, Han Y, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res 2021;49:D1420-30. [Crossref] [PubMed]
21. Breiman L. Random forests. Machine Learning 2001;45:5-32. [Crossref]
22. Ester M, Kriegel HP, Xu X. XGBoost: A scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. 2016:785-79. Geographical Analysis 2022.
23. Lundberg SM, Nair B, Vavilala MS, et al. Explainable machine-learning predictions for the prevention of hypoxaemia during surgery. Nat Biomed Eng 2018;2:749-60. [Crossref] [PubMed]
24. Xu L, Deng C, Pang B, et al. TIP: A Web Server for Resolving Tumor Immunophenotype Profiling. Cancer Res 2018;78:6575-80. [Crossref] [PubMed]
25. Subramanian A, Narayan R, Corsello SM, et al. A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles. Cell 2017;171:1437-1452.e17. [Crossref] [PubMed]
26. Võsa U, Claringbould A, Westra HJ, et al. Large-scale cis- and trans-eQTL analyses identify thousands of genetic loci and polygenic scores that regulate blood gene expression. Nat Genet 2021;53:1300-10. [Crossref] [PubMed]
27. Okła K, Farber DL, Zou W. Tissue-resident memory T cells in tumor immunity and immunotherapy. J Exp Med 2021;218:e20201605. [Crossref] [PubMed]
28. Mao X, Xu J, Wang W, et al. Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives. Mol Cancer 2021;20:131. [Crossref] [PubMed]
29. Sancar A, Van Gelder RN. Clocks, cancer, and chronochemotherapy. Science 2021;371:eabb0738. [Crossref] [PubMed]
30. Llovet JM, Castet F, Heikenwalder M, et al. Immunotherapies for hepatocellular carcinoma. Nat Rev Clin Oncol 2022;19:151-72. [Crossref] [PubMed]
31. Raoul JL, Forner A, Bolondi L, et al. Updated use of TACE for hepatocellular carcinoma treatment: How and when to use it based on clinical evidence. Cancer Treat Rev 2019;72:28-36. [Crossref] [PubMed]
32. Tong M, Che N, Zhou L, et al. Efficacy of annexin A3 blockade in sensitizing hepatocellular carcinoma to sorafenib and regorafenib. J Hepatol 2018;69:826-39. [Crossref] [PubMed]
33. Tanouti IA, Fellah H, El Fihry R, et al. Association of Peroxisome Proliferator-Activated Receptor Gamma Coactivator 1 Alpha Coding Variants with Hepatocellular Carcinoma Risk in the Moroccan Population: A Case-Control Study. Asian Pac J Cancer Prev 2023;24:3689-96. [Crossref] [PubMed]
34. Huang DQ, Mathurin P, Cortez-Pinto H, et al. Global epidemiology of alcohol-associated cirrhosis and HCC: trends, projections and risk factors. Nat Rev Gastroenterol Hepatol 2023;20:37-49. [Crossref] [PubMed]
35. Llovet JM, Montal R, Sia D, et al. Molecular therapies and precision medicine for hepatocellular carcinoma. Nat Rev Clin Oncol 2018;15:599-616. [Crossref] [PubMed]
36. Xuan W, Khan F, James CD, et al. Circadian regulation of cancer cell and tumor microenvironment crosstalk. Trends Cell Biol 2021;31:940-50. [Crossref] [PubMed]
37. Fagiani F, Di Marino D, Romagnoli A, et al. Molecular regulations of circadian rhythm and implications for physiology and diseases. Signal Transduct Target Ther 2022;7:41. [Crossref] [PubMed]
38. Reinke H, Asher G. Circadian Clock Control of Liver Metabolic Functions. Gastroenterology 2016;150:574-80. [Crossref] [PubMed]
39. Mukherji A, Bailey SM, Staels B, et al. The circadian clock and liver function in health and disease. J Hepatol 2019;71:200-11. [Crossref] [PubMed]
40. Charlot A, Hutt F, Sabatier E, et al. Beneficial Effects of Early Time-Restricted Feeding on Metabolic Diseases: Importance of Aligning Food Habits with the Circadian Clock. Nutrients 2021;13:1405. [Crossref] [PubMed]
41. Martin RA, Viggars MR, Esser KA. Metabolism and exercise: the skeletal muscle clock takes centre stage. Nat Rev Endocrinol 2023;19:272-84. [Crossref] [PubMed]
42. Chen W, Zheng D, Chen H, et al. Circadian Clock Regulation via Biomaterials for Nucleus Pulposus. Adv Mater 2023;35:e2301037. [Crossref] [PubMed]
43. Wang C, Zeng Q, Gül ZM, et al. Circadian tumor infiltration and function of CD8(+) T cells dictate immunotherapy efficacy. Cell 2024;187:2690-2702.e17. [Crossref] [PubMed]
44. Kinouchi K, Sassone-Corsi P. Metabolic rivalry: circadian homeostasis and tumorigenesis. Nat Rev Cancer 2020;20:645-61. [Crossref] [PubMed]
45. Saran AR, Dave S, Zarrinpar A. Circadian Rhythms in the Pathogenesis and Treatment of Fatty Liver Disease. Gastroenterology 2020;158:1948-1966.e1. [Crossref] [PubMed]
46. Jamshed H, Beyl RA, Della Manna DL, et al. Early Time-Restricted Feeding Improves 24-Hour Glucose Levels and Affects Markers of the Circadian Clock, Aging, and Autophagy in Humans. Nutrients 2019;11:1234. [Crossref] [PubMed]
47. Sun YF, Wu L, Liu SP, et al. Dissecting spatial heterogeneity and the immune-evasion mechanism of CTCs by single-cell RNA-seq in hepatocellular carcinoma. Nat Commun 2021;12:4091. [Crossref] [PubMed]
48. Costa A, Kieffer Y, Scholer-Dahirel A, et al. Fibroblast Heterogeneity and Immunosuppressive Environment in Human Breast Cancer. Cancer Cell 2018;33:463-479.e10. [Crossref] [PubMed]
49. Thorsson V, Gibbs DL, Brown SD, et al. The Immune Landscape of Cancer. Immunity 2018;48:812-830.e14. [Crossref] [PubMed]
50. Zha M, Tian T, Xu W, et al. The circadian clock gene Bmal1 facilitates cisplatin-induced renal injury and hepatization. Cell Death Dis 2020;11:446. [Crossref] [PubMed]
51. Petit E, Milano G, Lévi F, et al. Circadian rhythm-varying plasma concentration of 5-fluorouracil during a five-day continuous venous infusion at a constant rate in cancer patients. Cancer Res 1988;48:1676-9. [PubMed]
52. Shafi AA, McNair CM, McCann JJ, et al. The circadian cryptochrome, CRY1, is a pro-tumorigenic factor that rhythmically modulates DNA repair. Nat Commun 2021;12:401. [Crossref] [PubMed]
53. Tsang FH, Law CT, Tang TC, et al. Aberrant Super-Enhancer Landscape in Human Hepatocellular Carcinoma. Hepatology 2019;69:2502-17. [Crossref] [PubMed]
54. Miller S, Kesherwani M, Chan P, et al. CRY2 isoform selectivity of a circadian clock modulator with antiglioblastoma efficacy. Proc Natl Acad Sci U S A 2022;119:e2203936119. [Crossref] [PubMed]
55. van der Watt PJ, Roden LC, Davis KT, et al. Circadian Oscillations Persist in Cervical and Esophageal Cancer Cells Displaying Decreased Expression of Tumor-Suppressing Circadian Clock Genes. Mol Cancer Res 2020;18:1340-53. [Crossref] [PubMed]
56. Chan AB, Parico GCG, Fribourgh JL, et al. CRY2 missense mutations suppress P53 and enhance cell growth. Proc Natl Acad Sci U S A 2021;118:e2101416118. [Crossref] [PubMed]
57. Xue J, Yi J, Zhu X. Knockdown of UCHL3 inhibits esophageal squamous cell carcinoma progression by reducing CRY2 methylation. Hum Cell 2022;35:528-41. [Crossref] [PubMed]
58. Takaki T, Montagner M, Serres MP, et al. Actomyosin drives cancer cell nuclear dysmorphia and threatens genome stability. Nat Commun 2017;8:16013. [Crossref] [PubMed]
59. Chen KT, Pernelle K, Tsai YH, et al. Liver X receptor α (LXRα/NR1H3) regulates differentiation of hepatocyte-like cells via reciprocal regulation of HNF4α. J Hepatol 2014;61:1276-86. [Crossref] [PubMed]
60. Krauß L, Urban BC, Hastreiter S, et al. HDAC2 Facilitates Pancreatic Cancer Metastasis. Cancer Res 2022;82:695-707. [Crossref] [PubMed]
61. Hu Z, Dong L, Li S, et al. Splicing Regulator p54(nrb) /Non-POU Domain-Containing Octamer-Binding Protein Enhances Carcinogenesis Through Oncogenic Isoform Switch of MYC Box-Dependent Interacting Protein 1 in Hepatocellular Carcinoma. Hepatology 2020;72:548-68. [Crossref] [PubMed]
62. Liang Q, Wang Y, Lu Y, et al. RANK promotes colorectal cancer migration and invasion by activating the Ca(2+)-calcineurin/NFATC1-ACP5 axis. Cell Death Dis 2021;12:336. [Crossref] [PubMed]
63. Luo P, Feng X, Jing W, et al. Clinical and Diagnostic Significance of Homer1 in hepatitis B virus-induced Hepatocellular Carcinoma. J Cancer 2018;9:683-9. [Crossref] [PubMed]
64. Fang L, Lv J, Xuan Z, et al. Circular CPM promotes chemoresistance of gastric cancer via activating PRKAA2-mediated autophagy. Clin Transl Med 2022;12:e708. [Crossref] [PubMed]
65. LeBleu VS, O'Connell JT, Gonzalez Herrera KN, et al. PGC-1α mediates mitochondrial biogenesis and oxidative phosphorylation in cancer cells to promote metastasis. Nat Cell Biol 2014;16:992-15. Erratum in: Nat Cell Biol 2014;16:1125. [Crossref] [PubMed]