In this study, we analyzed two publicly available single-cell transcriptome datasets on liver cancer annotated with clinical information (18). We integrated and normalized the two datasets according to a previously reported pipeline (19). The merged single-cell transcriptome data were log-normalized and scaled for downstream analysis using Seurat V4.0.1. In addition to single-cell transcriptome data, we downloaded multiple bulk liver cancer transcriptome datasets and their corresponding clinical information to verify results obtained from the merged single-cell transcriptome datasets. The log2 normalized and transformed gene expression matrix of the Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) dataset was downloaded from the University of California, Santa Cruz (UCSC) Xena browser (https://xenabrowser.net/). Transcriptomic and clinical data from 159 pairs of tumor and normal samples taken from Chinese HCC patients with HBV infection were downloaded from the National Omics Data Encyclopedia database (https://www.biosino.org/node/) using the accession number OEP000321 (20). The LIRI-JP dataset, including transcriptome data from 231 HCC samples, was downloaded from the International Cancer Genome Consortium Data Portal (https://dcc.icgc.org/). Gene expression microarray data and detailed clinical information on GSE14520 (21) (including 221 HCC samples) were downloaded from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/).

CRD-related genes (CRGs) were collected from the CircaDB database (22). We employed gene set enrichment analysis to infer the CRD level of each cell based on the CRG expression matrix. First, we normalized all the sample matrices to minimize the effect of outliers. The relative ranking of CRDs among all genes reflects the overall CRD status of the sample. We used single-sample gene set enrichment analysis (ssGSEA) to calculate the relative CRD enrichment score of each sample (LiverCRD score), which represented the CRD status of the sample. We then employed a random sampling strategy to determine the abnormal gene set activity score. All genes were divided into 50 expression bins based on their average expression. The number of occurrences of CRGs, that is, their frequency within each bin, was calculated and identified as N. Based on random sample N timings, random signature genes from each bin were chosen. In other words, the total number of random signature genes K matched the number of CRGs. This process was repeated 1000 times at random. We then created a random score as the mean of the K x 1000 random signatures sampled to define the background level of CRGs. The threshold for determining abnormal CRD levels was set at 75% based on quartiles, which was -0.24.

We performed a weighted correlation network analysis (WGCNA) (23) to identify CRGs significantly associated with one or more clinical properties of the samples. This analysis was performed using the R package WGCNA V1.70-3. We screened the gene modules that were significantly positively correlated with the tumor stage (p< 0.05). Finally, 206 genes were screened for subsequent analysis and modelling.

We used data obtained from CellphoneDB (24) to analyze intercellular interactions. First, we investigated whether there was significant interaction (p< 0.05) between the two cell types (tumor and normal) by evaluating the average expression levels of annotated ligand-receptor pairs in the STRING database (25). Ligand-receptor pairs associated with p< 0.05 were subsequently screened to evaluate the relationship between the two cell types. Given that cytokines are critical for intercellular communication in cancer progression, we further analyzed cytokine signaling activity in the transcriptomic profile using CytoSig (26) to predict the effect of cytokines from the transcriptome-level signaling cascade.

The differentially expressed gene levels of samples in the high-CRD and low-CRD groups were calculated using the R package DESeq2 V1.30.1 (27). We set FDR< 0.05 and |log₂FC| > 1 as the criteria for differentially expressed genes. Pathway enrichment analysis based on cancer hallmarks and the Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling (28) pathway was conducted using the enricher function of the R package clusterProfile V3.18.1 (29).

We used CIBERSORT (30) to determine the relative composition of immune cells in TCGA-LIHC samples based on the gene expression matrix. The overall immune and stromal infiltration levels were estimated using the R package Estimate V1.0.13 (31).

We downloaded gene expression matrices and half-maximal inhibitory concentration (IC50) values of 805 cell lines treated with 198 drugs from the Genomics of Drug Sensitivity in Cancer (GDSC) database (32). Based on this training dataset, we used the R package oncoPredict V0.2 to build a ridge regression model and predict the response of each TCGA sample to each drug. The change in response of the high-risk and low-risk groups to each drug represented the difference in drug response between the two groups. Tumor immune infiltration data were downloaded from the Tumor Immune Dysfunction and Exclusion (TIDE) database (http://tide.dfci.harvard.edu/download/). Immunophenoscores (IPSs) were downloaded from the Cancer Immunome Atlas (https://tcia.at/home). Higher tumor-infiltrating cell exclusion scores and lower IPSs predicted poorer responses to immunotherapy.

Somatic mutation data on the TCGA-LUAD samples were downloaded from the UCSC Xena browser (https://xenabrowser.net/). The mutational landscape and somatic mutation interactions of samples with high and low CRD levels were analyzed using maftools V2.6.05. Among the various types of mutations, Missense, Nonsense, Frame_Shift_Ins, Frame_Shift_Del, In_Frame_Ins, In_Frame_Del, Splice_Site, Translation_Start_Site, and Nonstop were considered nonsynonymous mutant variants. Silent and other types of mutations, including introns, 3’ untranslated regions, 5’ untranslated regions, and intergenic regions, were considered synonymous mutant variants. Synonymous mutations were also considered wild-types because they did not affect proteins. Frameshift and nonsynonymous mutant variants were encoded as truncating mutations, consistent with previous studies (33). Mutation proportions were assessed using one-sided Z-tests and two-sided χ² tests. Statistical significance was set at p< 0.05.

To better apply CRD scores in clinical sample analysis and make it possible to evaluate patient prognosis, we used the LASSO-Cox regression model to screen potential genes significantly associated with prognosis based on the 206 CRGs that were significantly positively correlated with the tumor stage. After 1,000 permutations and cross-validation, 14 genes were screened. Linear combinations of these genes (based on their expression levels) were used to calculate the CRD scores of the HCC samples from each patient. The minimum criterion determined the regression coefficient. The final CRD risk score of each sample was calculated as follows:

We subsequently divided the samples into high-risk and low-risk groups based on the median CRD score. The statistical significance of overall survival (OS) was determined using the log-rank test and visualized as Kaplan–Meier curves using the R package survminer V0.4.9.

Most patients in both single-cell transcriptome datasets received PD-1 or PD-L1/CTLA-4 monoclonal antibody therapy and predominantly exhibited stage IV liver cancer development ( Figure 1A ). The cell numbers identified for each patient ranges from 109 to 1, 416. The merged, regrouped, and annotated single-cell transcriptome datasets are shown in Figures S1A, B . As shown in Figures 1B, C , malignant cells differed significantly between patients, whereas immune cells were more consistent. The epithelial score of malignant cells was markedly higher than that of stromal and immune cells ( Figure S1C ), which is consistent with previous studies (34). Marker genes in different clusters were highly expressed in specific immune cells ( Figure S1D ), further supporting the accuracy of our clustering. Based on the collected CRGs ( Table S1 ), the CRD score of malignant cells was significantly higher than that of other cells, indicating that the malignant cell cycle was more severely disordered ( Figure 1D ). We further divided the malignant cells into groups with high, medium, and low CRD scores ( Figure S1E ). We found that cells with a high CRD score were enriched in various metabolic signaling pathways. In contrast, cells with a low CRD score were significantly enriched in signaling pathways such as cell junctions ( Figure 1E ). Moreover, the CRD score of malignant cells varied widely between different samples ( Figure 1F ), indicating that liver cancer has high intratumoral heterogeneity.

We subsequently screened for CRGs associated with clinical factors. First, we used the WGCNA method to analyze the single-cell transcriptome matrix ( Figure S2A ) and divide the CRGs into multiple modules ( Figure 1G ). The relationship between core genes and clinical variables in each module was then calculated ( Figure 1H ), and genes significantly associated with these variables were identified ( Table S2 ). We found that these genes were greatly enriched in tumorigenesis and development-related signaling pathways, such as epithelial-mesenchymal transition, MAPK signaling, P53 signaling, hypoxia, and cytokine receptor interaction ( Figure S2B ), suggesting that these CRD scores are closely related to the malignant progression of HCC. In addition, the CRD score was significantly correlated with clinical variables, such as age, sex, cancer stage, and the viral infection status of the samples ( Figure 1I ).

Cytokine interactions between immune cells and malignant cells in the tumor microenvironment (TME) have always been the focus of research in antitumor immunotherapy (35, 36). Owing to the significant enrichment of cytokine interaction signaling pathways reported in Section 3.1, we speculated that cells in the high-CRD and low-CRD groups and immune cells in the TME are also involved in other types of intercellular communication. We thus used the CellphoneDB to analyze intercellular interactions and construct an intercellular communication network ( Figure 2A ). Notably, among the chemokines, high-CRD and low-CRD groups exhibited distinct intercellular communication with immune cells ( Figure 2B ), including CCL20, CXCL12, and CCL5, all of which were affected by CRD. In the low-CRD group, the accumulation of CCL5/CCR5, CCL5/ACKR1, and CXCR3/CCL20 suggests the accumulation of CD8(+) T cells. These results indicate that the low-CRD group may have had a higher level of immune infiltration. In comparison, the high-CRD group may have had a lower level of immune infiltration and inconsistent immunotherapy responses. We also compared immune co-suppressive factors ( Figure 2C ) and immune co-stimulatory factors ( Figure 2D ). However, we did not find significant differences between the low-CRD and high-CRD groups. To our knowledge, in the tumor microenvironment, chemokines can be expressed in various types of cells, and their expression levels are high, while stimulatory factors and inhibitory factors are not much detectable, and their expression levels are also limited. Therefore, chemokine expression is more easily detected and differential chemokines are more easily identified, so it is easier for us to find differences in chemokines from the results.

Cytokines are critical for intercellular communication in cancer progression (37). To investigated whether CRDs regulate cytokine signaling activity, we calculated the cytokine activities at the single-cell level. As shown in Figure 2E , cytokines such as IL17A and IL1A were upregulated in TAM, and cytokines such as MCSF, IL3, IGF1, and INS were significantly upregulated in the high-CRD group. IL17A and IL1A are two well-known regulators of macrophages (38, 39), indicating the correctness of the cytokine activity we calculated. MCSF is the macrophage colony-stimulating factor that regulate the differentiation of myeloid lineage cells (40), and IL3 plays functional role on regulating the myofibroblastic differentiation (41). The diversity of CRD scores suggests differences in the TME composition—CRD may be a significant regulator of immune cell accumulation in HCC.

Based on the CRGs related to clinical variables screened in Section 3.1, we used the ssGSEA method to evaluate the CRD status of the bulk samples. We defined the enrichment score as the LiverCRD score ( Table S3 ), which significantly correlated with tumor stage ( Figures 3A , S3A ). We further divided the samples into high and low median scores and compared their gene expression and pathway activity levels ( Figure 3B and Table S4 ). Similar to the single-cell transcriptome results in Section 3.1, in both the cancer hallmark and KEGG signaling pathway datasets, we observed the differential activation of many metabolism-related signaling pathways, such as bile acid and fatty acid metabolism. This coincided with the activation of many signals related to tumor development, such as the MYC signaling pathway, E2F target, oxidative phosphorylation, interferon response, and mTOR signaling pathway ( Figures 3C, D ).

These results suggest that CRD is likely associated with the remodeling of the immune microenvironment. Although the LiverCRD and immune/stromal infiltration scores estimated using the Estimate algorithm did not significantly correlate with each other ( Figure S3B ), the two groups exhibited significant differences in T cell and macrophage levels ( Figure 3E ). Many immune checkpoints were higher in the high-LiverCRD group ( Figure S4A ), and the expression of immune checkpoints, such as PD-1 and CTLA-4, was significantly positively correlated with the LiverCRD score ( Figure S4B). Many immune-related molecules, such as MHC molecules, immune stimulators, and immune inhibitors, were regulated by CRD ( Figure 3F ). Based on immune cell exclusion and dysfunction levels downloaded from the TIDE database, the high-LiverCRD group had a higher level of immune cell rejection ( Figure S4C ). Our predictions suggested that the high-LiverCRD group had worse response to immunotherapy than the low-LiverCRD group ( Figures 3G , S4D ). These results reflect those of a previous single-cell transcriptome analysis, which showed that the immune microenvironment of the high-CRD group had a lower level of immune infiltration and lower CD8(+) T cell content than the low-CRD group (42). Based on the gene expression matrix and IC50 value of the GDSC database, we constructed a ridge regression model to predict the drug response of TCGA samples ( Table S5 ). The high-LiverCRD group responded better to MDM2, MEK, and mTOR inhibitors than other drugs ( Figure 3H ), reflecting the higher activity of these pathways in the high-CRD group in the previous analysis.

The precise regulation of the cell cycle relies on a series of regulatory molecules, among which mutations are particularly critical (43). We thus analyzed the mutational landscape of the high-CRD and low-CRD groups ( Figure 4A ). As shown in Figure 4B , TP53 and MYO18B more frequently mutated in the high-CRD group than the low-CRD group. In addition, TP53 mutations in the high-CRD group were more likely to be enriched via truncating mutations that affect protein function. In contrast, the number of TP53 mutations in the low-CRD group that did not affect protein function was relatively high ( Figure 4C ), suggesting that CRD may be regulated by TP53 mutation status. Notably, although ARID1A mutation levels and types were not significantly different between the two groups ( Figure 4D ), ARID1A mutations were associated with lower patient survival rates in the high-CRD group than in the low-CRD group ( Figure 4E ), suggesting that ARID1A mutations in the former were affected by factors performing different functions. We also compared the somatic mutation interactions between the two groups. Although the two groups of mutations generally co-occurred, the specific crosstalk differed between them ( Figure 4F ). The ARID1A/VCAN, KMT2C/UBR4, and LAMA1/CSMD1 pair mutations co-occurred in the high-CRD group. The MUC2/AHNAK, PRKDC/SLIT2, and MUC16/MUC2 pairs of mutations co-occurred in the low-CRD group.

After 1,000 permutations and cross-validation, we determined the parameters of the LASSO-Cox regression model ( Figure 5A ). The model was constructed using 14 genes and calculated their respective regression coefficients ( Figure 5B and Table S6 ). As shown in Figure 5C , patients who died had significantly higher CRD risk scores, indicating that CRD risk significantly affected patient survival ( Figure 5D ). As shown in Figure 5E , the median survival times of the high-risk and low-risk groups were approximately 800 and 2,500 days, respectively. We verified the robustness of the high-risk scores using three independent datasets. As shown in Figure 5F , the CRD risk scores of different datasets significantly correlated with the hazard ratio, and their respective patient survival curves also differed significantly ( Figure S5 ). Although these datasets were sequenced by various countries, the survival analysis was performed within a single dataset, and survival differences were consistent across datasets. These results suggest that the CRD risk score is a robust indicator of survival in HCC patients.

Compared with clinical factors, the CRD score could assess patient risk similarly to tumor stage ( Figure 5G ). Time-dependent receiver operating characteristic curve analysis also indicated that the CRD risk score was effective in assessing patients 1-, 3-, and the effect was stable in the 5-year survival state ( Figure 5H ). In addition, we constructed a nomogram for multivariate analysis, showing that the CRD risk score predicted overall survival ( Figure S6A ). Calibration curves showed that the CRD risk score did not deviate significantly from the actual value in predicting one-, three-, and five-year patient survival ( Figures S6B–D , respectively). These results suggest that the CRD risk score is a suitable indicator for monitoring HCC patient survival.