Introduction

In 2020, the global incidence of primary liver cancer reached 906,000 new cases, ranking it as the sixth most prevalent malignancy. Despite being the sixth most common form of cancer worldwide, liver cancer ranks as the third highest cause of cancer-related deaths.1 Hepatocellular carcinoma (HCC) constitutes around 80–90% of all primary malignant liver tumors.2 Recognized risk factors for the development of HCC include infection with hepatitis B virus, infection with hepatitis C virus, liver cirrhosis, excessive alcohol consumption, non-alcoholic fatty liver disease, type II diabetes, and obesity.3,4 HCC presents a notably severe scenario in China, ranking fourth in malignant tumor incidence and second in mortality.5,6 Specifically in regions like Guangxi, prolonged exposure to hepatitis B and aflatoxins has led to the sustained dominance of liver cancer in both incidence and mortality among all malignancies in China.7 The main treatment for HCC is surgical resection. However, the 5-year survival rate of postoperative patients is only 50% to 68%.8 Recently, the introduction of tyrosine kinase inhibitors (TKIs) such as Sorafenib and Lenvatinib, along with immune checkpoint inhibitors (ICIs) like Pembrolizumab and Atezolizumab, has expanded the therapeutic choices for patients with HCC. Nevertheless, the effectiveness of TKI or ICI monotherapy for HCC alone is considered unsatisfactory. The combination of ICIs and TKIs can elevate the objective response rate to 30%-40%.9,10 When further combined with local treatments, such as transarterial chemoembolization (TACE) or hepatic arterial infusion chemotherapy, the objective response rate can exceed 60%.11 Various biomarkers, including PD-L1 expression levels, nutritional status, treatment-related adverse events, neoantigens, gut microbiota, tumor mutational burden (TMB), and immune cell status, can to some extent predict the efficacy of ICI therapy.12–16 The immune status of cancer can function as a predictor for both tumor prognosis and treatment.

The tumor microenvironment is an intricate amalgamation of cellular components, primarily comprising immune cells, inflammatory cells, cancer-associated fibroblasts, extracellular matrix, and the network of blood and lymphatic vessels. Notably, immune cells play a crucial role in shaping the response to targeted and immunotherapy.17 Dendritic cells (DCs) are crucial immune-regulating cells with the ability for antigen presentation, serving as sentinel cells in the immune system and capable of activating T lymphocytes.18 DCs can also contribute to immune tolerance.19 Regulatory DCs, constituting approximately 15% of circulating monocytes in HCC patients, promote systemic immune suppression by producing high levels of IL-10.20 Recently, there has been an increasing recognition that the enrichment and functionality of DCs can influence the efficacy of diverse cancer therapies. Chemotherapy-induced immunogenic cell death, under appropriate conditions, has the potential to amplify immunogenicity. For instance, the translocation of intracellular calcium-binding proteins to the surface of dying tumor cells can enhance the phagocytic activity of CD11c+ DCs.21 Research indicates that tumor cells treated with TKIs can induce the complete maturation of DCs.22 TKIs could also inhibit the ERK/β-catenin pathway, thereby promoting the recruitment and differentiation of DCs.23 PD-L1 on DCs has been recognized as a direct target for immunotherapy and a determinant of ICIs responses. By blocking PD-L1 on the surface of DCs, the interaction between CD80 and CD28 is facilitated, promoting the activation of T cells.24

This study assesses DCs infiltration across diverse HCC databases and constructs a dendritic cell-related signature (DCRS) model employing Cox regression and least absolute shrinkage and selection operator (LASSO) regression analyses. These genetic characteristics may have the capacity to distinguish between specific HCC subtypes and predict extended prognosis, immune cell infiltration levels, drug sensitivity, and treatment responsiveness of HCC patients.

Materials and Methods Public Datasets Collection

RNA sequencing data, gene mutation data, and clinical prognosis information for HCC patients were obtained from The Cancer Genome Atlas (TCGA) database, comprising 374 HCC tissues and 50 adjacent non-cancerous tissues (Table S1). Simultaneously, transcriptome data and clinical-pathological information for HCC from the GSE76427 dataset (115 samples) and GSE14520 dataset (242 samples) were acquired from the Gene Expression Omnibus database. Batch correction was performed using the ComBat algorithm provided by the sva R package across the TCGA, GSE76427, and GSE14520 datasets. Furthermore, datasets associated with HCC treated with TACE (GSE104580, 147 samples), Sorafenib (GSE109211, 140 samples), and anti-PD-1 immunotherapy (GSE202069, 17 samples), along with a single-cell dataset for HCC (GSE140228), were downloaded.

Clinical Samples Collection

With explicit written consent, we acquired 13 pairs of HCC tissues and paracancer tissues after hepatectomy between 2020 to 2023, previously received to anti-PD-1 (Camrelizumab) and TKI (Apatinib) therapies, from the Department of Hepatobiliary Surgery at the First Affiliated Hospital of Guangxi Medical University, named “Guangxi cohort”. In liver tissues, we isolated both tumor and paracancer tissues, each approximately the size of a soybean. The tissues were then preserved in an RNA stabilization solution and stored at −80°C. The study received ethical approval, and the assigned approval number is: [2024-E127-01].

DC-Related Gene and Consensus Clustering Analysis

Utilizing the TIME method, we evaluated immune infiltration in HCC samples from TCGA-LIHC, GSE76427, and GSE14520 datasets. The Weighted Gene Co-expression Network Analysis (WGCNA) was employed to explore co-expressed genes of DCs across these datasets. Venn diagrams were utilized to identify gene intersections, forming a DC-related gene set. Unsupervised hierarchical clustering was conducted on the DC-related gene set in the TCGA-LIHC dataset using the “ConsensusClusterPlus” R package. Differential expression genes (DEGs) analysis was performed using the limma package in R. Genes with a Benjamini–Hochberg-corrected p-value (adjusted p value) < 0.05 and |log2 fold change| > 1 were considered differentially expressed.

Functional Enrichment Analysis

Gene set enrichment analysis (GSEA), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) were employed to assess the biological roles of candidates using the “clusterProfiler” R package.

Identification of a DCRS for HCC

Based on unsupervised clustering results, the TCGA-LIHC dataset was categorized into two subtypes: cluster 1 and cluster 2. A candidate gene set was obtained by intersecting DEGs between cluster 1 and cluster 2 and those between HCC and adjacent tissues in TCGA-LIHC. The caret package was used to randomly split TCGA-LIHC 1:1 into training and validation sets. In the training set, univariate Cox regression was initially performed, followed by LASSO regression analysis using the glmnet package. The resulting genes with the minimum λ value were used to construct a multivariate Cox regression model with the smallest Akaike information criterion (AIC) through stepwise regression (backward). The risk score was calculated as follows: where denoted the expression level of the i-th dendritic cell-related gene, and represented its corresponding regression coefficient. The summation was over all n genes included in the risk score model. The gene set was categorized into high and low-risk groups based on the median risk score. The same procedures were applied to validate the risk score model in TCGA validation set, GSE76427, and GSE14520 datasets.

Construction of Clinical Prognosis Nomogram Model

To explore the relationship between the risk scores and clinical-pathological parameters, we examined the differences in risk scores among different clinical parameter groups. Univariate and multivariate Cox regression analyses were conducted in the TCGA dataset to identify independent prognostic parameters. A nomogram was constructed using the rms package, and its predictive accuracy was evaluated using calibration curves. Decision curve analysis was performed to determine the clinical application value of the prediction model, based on the ggDCA package.

Genomic Mutation Analysis

Utilizing TCGA mutation data, the Maftools package was employed to visualize mutation profiles for high and low risk groups. The TMB for each HCC patient was also calculated.

Immune Cell Infiltration

Three algorithms from the R package IOBR, TIMER, ESTIMATE, and xCell were used to determine the proportion of immune infiltrating cells based on the TCGA-LIHC dataset. TIMER uses an inverse convolution approach for estimating the proportion of six immune cell types. The xCell conducts cell type enrichment analysis using data on the 64 immune and stromal cell types’ gene expression. The ESTIMATE algorithm was employed to determine the immune score, tumor purity, stromal score and ESTIMATE score for tumors. Gene set variation analysis was conducted using the tumor microenvironment-related gene sets integrated in the R package IOBR.

Drug Sensitivity Prediction

Drug sensitivity analysis was conducted using the oncoPredict R package. Expression data of model genes were used to predict the sensitivity (IC50 values) of 367 drugs in the GDSC1 database. The differences in IC50 values between risk groups were compared, and the correlation between IC50 and risk score was assessed using Spearman correlation coefficient.

Guangxi Cohort RNA-Sequencing

RNA sequencing was performed on cancer and paracancerous tissues from patients undergoing hepatectomy. RNA was extracted from tumor and corresponding paracancerous samples utilizing the RNeasy Mini kit (QIAGEN). After quantitative analysis, rRNA was eliminated using the NEB Next rRNA Depletion Kit (NEB, E6310S), followed by reverse transcription to produce cDNA. Fastqc (v0.11.5) software was employed to evaluate raw sequencing data, and AdapterRemoval (ver.2.2.2) software was utilized to obtain clean reads by removing adapters and low-quality sequences. Clean reads were aligned to the hg19 complete transcriptome sequences through STAR (2.5.3a) software, and RSEM (1.3.0) software was employed for gene expression quantification.

Predicting the Response to Therapy

To assess the predictive potential of DCRS in HCC treatment, we analyzed the differences in risk scores between treatment response groups in GSE104580, GSE109211, and GSE202069 datasets using Wilcoxon test. Additionally, we validated the predictive ability of the risk score model using transcriptome data from the immunotherapy and targeted therapy-treated Guangxi cohort.

Single Cell RNA Sequencing Data Analysis

From the GSE140228 dataset, HCC and normal liver tissue samples were extracted. Harmony was used to remove batch effects, and SingleR was employed for cell annotation. Cell clustering was performed using the UMAP method based on the Seurat package.

Statistical Analysis

All statistical analyses were performed using R software (v4.3.1). The Shapiro–Wilk test was used to assess the distribution of continuous data. Data with a normal distribution were presented as mean ± standard deviation and analyzed using the t-test. Continuous data with a non-normal distribution were expressed as median (P25, P75) and compared using the Wilcoxon test and Kruskal–Wallis test. Categorical data were reported as frequencies (percentages) and analyzed with the chi-square test or Fisher’s exact test. Spearman’s rank correlation coefficient was used to evaluate relationships between continuous variables. Kaplan-Meier survival analysis and Log rank tests were conducted to assess survival differences, and the time-dependent receiver operating characteristic (ROC) curve was used to evaluate predictive performance. Statistical significance was defined as a P-value < 0.05. The Benjamini-Hochberg method was applied to correct for false discovery rates.

Results Identification of a DC-Related Gene Set

Based on WGCNA, in the TCGA cohort, the green module (919 genes) exhibited the highest correlation with DCs infiltration (cor = 0.99, P < 1e-200) (Figure 1A and D, Table S2). In GSE14520, the magenta module (550 genes) showed the highest correlation with DCs (cor = 0.96, P < 1e-200) (Figure 1B and E). In GSE76427, the blue module (1423 genes) demonstrated the highest correlation with DCs (cor = 0.95, P < 1e-200) (Figure 1C and F). The intersection of these three modules was obtained based on the venn diagram, resulting in 320 DC-related genes (Figure 1G). GO enrichment analysis and KEGG enrichment analysis were conducted on the DC-related gene sets, revealing enrichment in pathways such as leukocyte cell-cell adhesion, leukocyte-mediated immunity, and regulation of T cell activation, aligning with the functional aspects of DC antigen presentation (Figure 1H and I). The protein-protein interaction results indicated that CD4, PTPRC, and IL1, among others, serve as hub genes in the DC-related gene set (Figure S1).

Figure 1 Identify DC-related genes. (A-C), Use weighted gene co-expression network analysis to identify gene modules associated with immune cell infiltration in the TCGA (A), GSE14520 (B), and GSE76427 (C) datasets separately. The “r” represents the correlation coefficient; (D-F), Identify gene modules most correlated with DCs in the TCGA (D), GSE14520 (E), and GSE76427 (F) datasets; G, Venn diagrams to obtain DC-related genes based on the intersection. (H-I), Conduct Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses of DC-related genes.

Abbreviations: DC, dendritic cell; TCGA, The Cancer Genome Atlas.

Identification of DC-Related Clusters

Unsupervised clustering analysis was performed on all the 320 DC-related genes, resulting in the segregation of the TCGA dataset into two DC-related clusters (Figure 2A and B). Principal component analysis disclosed clear separation trends between the two clusters (Figure 2C). The heatmap illustrated a generally higher expression level of DC-related genes in cluster 1 compared to cluster 2 (Figure 2D). Additionally, we noted a more favorable overall survival (OS) in HCC patients from cluster 1 compared to cluster 2 (P = 0.037) (Figure 2E). Subsequent analysis of DEGs between cluster 1 and cluster 2 identified 1249 genes with significant differences (adjusted P < 0.05, |log2 fold change| > 1) (Figure 2F). GSEA was conducted on the DEGs in cluster 1 and cluster 2 (Figure 2G and H). The pathway curve plots illustrating the top 10 enriched pathways revealed enrichment in immune-related functional pathways such as CTL pathway, CTLA4 pathway, and immunotherapy by PD-1 blockade, indicating differences in the immune function between the two DC-related clusters.

Figure 2 Construct and analyze DC-related clusters. (A-B), Unsupervised clustering analysis identifies two DC-related clusters in the TCGA dataset; (C), Principal component analysis (PCA) demonstrates the separation trend of the two DC-related clusters; (D), Heatmap shows the enrichment levels of DC-related genes in the two clusters; (E), The Kaplan–Meier curve between the two DC-related clusters. (F), Volcano plot of differentially expressed genes between the two DC-related clusters; (G-H), Gene set enrichment analysis of differentially expressed genes.

Abbreviations: DC, dendritic cell; C1, cluster 1; C2, cluster 2; OS, overall survival, TCGA, The Cancer Genome Atlas.

Building and Validating DC-Related Signature

To identify characteristic genes of DC-related clusters, we intersected 1249 significantly DEGs with those differing between cancer and adjacent tissues in the TCGA dataset, resulting in 808 candidate genes (Figure 3A). Subsequently, patients in the TCGA database were randomly divided into a 1:1 ratio training set and a validation set. In the TCGA training set, Cox regression analysis on the 808 candidate genes identified 37 genes significantly associated with OS (P < 0.05) (Figure S2A). These 37 genes underwent LASSO regression and stepwise regression, resulting in a model with the minimum AIC (Figure 3B and C). Finally, 13 genes were included in the DCRS model. The risk score was calculated using the formula: Risk score = CLEC1B × −0.203 + CRIP1 × −0.304 + LTF × −0.335 + SLC4A10 × −0.226 + ALLC × −0.317 + PLA2G7 × 0.206 + CYP19A1 × 0.103 + CHGA × 0.195 + NR0B1 × 0.143 + SLC6A20 × 0.203 + CHRDL1 × 0.239 + CYP17A1 × −0.142 + SLC34A2 × 0.103 (Figure S2B). Additionally, Kaplan-Meier curve analysis was employed to evaluate the impact of high and low expression of DCRG genes on OS (Figure S3). Using the median of the risk score, the dataset was divided into high-risk and low-risk groups. In the TCGA training set, the low-risk group showed significantly better OS (P < 0.0001) (Figure 3D). Similarly, in the TCGA validation set (P = 0.0026), GSE76427 (P = 0.0051), and GSE14520 (P = 0.0047), the OS of the low-risk group surpassed that of the high-risk group (Figure 3E–G). The risk curves illustrated higher survival rates in the low-risk group (Figure 3H–K). ROC curves indicated that the area under curve (AUC) values of the risk scores in the TCGA training set at 1 year, 3 years, and 5 years were 0.873, 0.852, and 0.815, respectively (Figure 3L). Despite a slight decrease in the AUC values of the three validation sets, they still indicated the predictive potential of the risk score for the short-term and long-term survival of HCC patients (Figure 3M–O).

Figure 3 Constructing and validating the DC-related signature. (A), Venn diagram identifies the intersection of differentially expressed genes in DC-related clusters and differentially expressed genes between TCGA tumor and adjacent non-tumor tissues; (B-C), Construction of the DC-related signature based on LASSO regression; (D-G), Survival curves of DCRS risk groups in the TCGA training set, TCGA validation set, GSE76427, and GSE14520; (H-K), Risk score scatter plots of the TCGA training set, TCGA validation set, GSE76427, and GSE14520; (L-O), Time‐dependent receiver operating characteristic curves of the risk scores in the TCGA training set, TCGA validation set, GSE76427, and GSE14520.

Abbreviations: DC, dendritic cell; C1, cluster 1; C2, cluster 2; DEG, differentially expressed gene; OS, overall survival; AUC, area under curve; TCGA, The Cancer Genome Atlas.

Association of DCRS with Clinical Parameters and Nomogram

In the TCGA dataset, the risk scores were significantly lower in the T1 stage compared to both T2 (P=0.0002) and T3 stages (P=0.0008), as well as T4 (P=0.007). Moreover, in stage I, HCC patients had lower risk scores compared to patients in stages II (P=0.0005) and III (P=0.0002) (Figure 4A and B). Pathological grade reveals that the risk scores in grade 1 was significantly lower than in grade 2 (P=0.0022), grade 3 (P=2.2e-07) and grade 4 (P=0.0059) of HCC (Figure 4C). The risk scores exhibited higher expression in advanced tumors than in early tumors, consistent with significantly elevated scores in deceased patients compared to surviving patients (P=1.9e-09) (Figure 4D). When assessing the predictive capability of the risk score using ROC curves, the AUC values were 0.873 in TCGA training set, 0.725 in TCGA validation set, 0.626 in GSE76427, and 0.673 in GSE14520, slightly higher than other clinical parameters (Figure 4E–H). Based on univariate and multivariate Cox regression analysis in the TCGA dataset, the DCRS risk group and tumor staging acted as independent prognostic factors (Figure 4I and J). Constructing a model with parameters such as risk group, stage, pathological grade, age, and gender, the nomogram suggested that the risk group contributes the most substantially to predicting prognosis (Figure 4K). Calibration curves exhibited good consistency between predicted and actual outcomes at 1, 3, and 5 years in the nomogram (Figure 4L). Clinical decision curves indicated that the nomogram model provides greater benefits compared to other clinical parameters (Figure 4M).

Figure 4 Dendritic cell-related signature association with clinical parameters and construction of a nomogram. (A-D), Boxplots of risk scores for hepatocellular carcinoma patients stratified by different T stages, tumor stages, pathological grades, and survival status. (E-H), Receiver operating characteristic curves of risk scores and clinical parameters in the TCGA training set, TCGA validation set, GSE76427, and GSE14520 datasets. (I-J), Univariate and multivariate Cox regression analyses in the TCGA-LIHC dataset. (K), Nomogram combining risk groups with clinical parameters. (L), Calibration curves of nomogram model at 1, 3, and 5 years in the nomogram. (M), Decision curves of the nomogram and clinical parameters.\.

Abbreviations: AUC; area under the curve; TNM, Tumor Node Metastasis; HR, hazard ratio; OS, overall survival; TCGA, The Cancer Genome Atlas.

Association of DCRS with Mutational Landscape and TMB Analysis

The distinct molecular subtypes of tumors reflect their heterogeneity. Therefore, we further explored the mutational features of high and low risk groups. Using TCGA mutation data, we delineated the top 20 genes in the mutational landscape of both groups. Overall mutation frequencies were comparable between the two groups (86.63% vs 85.63%). Notably, TP53 mutations (40%) dominated in the high-risk group, while CTNNB1 (29%) mutations prevailed in the low-risk group (Figure 5A and B). Moreover, the risk scores in TP53-mutated HCC patients significantly exceeded that in TP53 wild-type patients (P=3.9e-08), whereas the risk scores in CTNNB1-mutated HCC patients were lower than in CTNNB1 wild-type patients (Figure 5C and D). Subsequently, TMB was calculated for the TCGA population, revealing a significantly higher TMB in group cluster 2 than in cluster 1 (P = 0.0005). However, no statistically significant difference in TMB was observed between high and low-risk groups (Figure 5E and F). Joint survival analysis of DCRG risk and TMB groups indicated a better prognosis in the low-risk group within the TMB subgroup (P <0.0001) (Figure 5G). Nevertheless, no statistically significant difference was found in the joint survival analysis between DC-related clusters and TMB groups (Figure 5H).

Figure 5 Analysis of gene mutations and TMB in DCRS risk groups. (A-B), Mutation landscape of top 20 genes in DCRS high-risk and low-risk groups. (C), Boxplot of risk scores for TP53 wild-type and mutation types. (D), Boxplot of risk scores for CTNNB1 wild-type and mutation types. (E), Boxplot of TMB in DC-related clusters. (F), Boxplot of TMB in DCRS risk groups. (G), Joint survival curve of risk groups and TMB. (H), Joint survival curve of DC-related clusters and TMB.

Abbreviations: TMB, tumor mutational burden. DCRS, dendritic cell-related signature; C1, cluster 1; C2, cluster 2; OS, overall survival.

Immune Infiltration Analysis

Gene set variation analysis reveals enrichment of gene sets related to tumor-associated macrophages (TAM), T cell exhaustion, T cell regulatory, and ICB resistance in the low-risk group, suggesting potential low sensitivity to ICIs (Figure 6A and B). Evaluation of tumor immune infiltration using the xCell algorithm demonstrated higher enrichment scores for CD8+ T cells, DCs, PDCs, ADCs, and CDCs in C1 compared to C2 (Figure 6C). Additionally, the low-risk group exhibited higher enrichment scores for CD8+ T cells, PDCs, ADCs, and CDCs compared to the high-risk group (Figure 6D). ESTIMATE algorithm assessment revealed that the stromal score in the low-risk group was significantly higher than that in the high-risk group (Figure 6E). Furthermore, both stromal score and immune score in cluster 1 were significantly higher than in cluster 2 (Figure 6F). Apart from CYP17A1, the expression of other DCRS genes is positively correlated with the infiltration of B cells, DCs, macrophages, neutrophils, CD4 T cells, and CD8 T cells.(Figure 6G). The IPS in the low-risk group was significantly higher than in the high-risk group, indicating a strong immunogenicity in the low-risk group (Figure S4).

Figure 6 Immune infiltration analysis. (A-B), Gene set variation analysis of DC-related clusters and risk groups; (C-D), Immune infiltration analysis based on the xCell algorithm in DC-related clusters and risk groups. The red text on the x-axis indicated dendritic cells and their subpopulations.; (E-F), Immune infiltration analysis based on the ESTIMATE algorithm in DC-related clusters and risk groups; G, Heatmap showing the correlation between DC-related signature genes and immune infiltrating cells.

Abbreviations: DC, dendritic cell; GPAGs: good prognosis angiogenesis genes; PPAGs, poor-prognosis angiogenesis genes; Pan-F TBRs, panfibroblast TGFβ response characteristics; MDSC, myeloid-derived suppressor cells; CAF, cancer associated fibroblast; EMT, epithelial-mesenchymal transition; ICB, immune checkpoint blockade; TME, tumor microenvironment; GEP, gene expression profile; TCR, T cell receptor; BCR, B cell receptor; CMP, common myeloid progenitor; MPP, multipotent progenitor; NK, natural killer cell; Tem, effector memory T cell; Tgd cells, T gamma delta cells; GMP, granulocyte-monocyte progenitor; Tregs, regulatory T cells; Tcm, central memory T cell; pDC, plasmacytoid dendritic cell; NKT, natural killer T cell; CLP, common lymphoid progenitor; MSC, mesenchymal stem cell; aDC, activated dendritic cell; Th1, T helper 1 cell; Th2, T helper 2 cell; iDC, immature dendritic cell; MEP, megakaryocyte-erythrocyte progenitor; cDC, conventional dendritic cell; HSC, hematopoietic stem cell.

Notes: *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001.

Predictive Potential of DCRS for HCC Treatment Response

To assess the potential of DCRS in predicting the efficacy of HCC treatment, we compared DCRS risk scores across different HCC treatment datasets. In the TACE-treated HCC dataset GSE104580, the risk scores in the response group were lower than in the non-response group (P = 0.01), consistent with the bar chart suggesting a higher response rate in the low-risk group (Figure 7A). In the Sorafenib-treated HCC dataset GSE109211, the risk scores in the response group were higher than in the non- response group (P = 0.041) (Figure 7B). In the PD1-treated HCC dataset GSE202069, the risk scores in the response group was higher than in the non-response group (P = 0.027) (Figure 7C).

Figure 7 Therapeutic response prediction and drug sensitivity analysis. (A), Boxplots and bar graphs of risk scores in the GSE104580 dataset treated with TACE; (B), Boxplots and bar graphs of risk scores in the GSE109211 dataset treated with Sorafenib; (C), Boxplots and bar graphs of risk scores in the GSE202069 dataset treated with anti-PD1 therapy; (D-O), Boxplots and correlation scatter plots of drug sensitivity analysis for different drugs.

Abbreviation: TACE, transarterial chemoembolization.

Drug Sensitivity Prediction

To further investigate the potential application of DCRS in personalized HCC treatment, we evaluated the IC50 values of different drugs between high and low risk groups. Drug sensitivity results indicated that the IC50 values of imatinib, cetuximab, paclitaxel, vinorelbine, gemcitabine, and epothilone were higher in the low-risk group, suggesting that high-risk patients may benefit more from these drugs (Figure 7D–I). In contrast, the IC50 values of doxorubicin, phenformin, tretinoin, capivasertib, dactolisib, and olaparib were lower in low-risk group, indicating that low-risk patients may benefit more from these therapies (Figure 7J–O). The drug sensitivity prediction of cluster 1 and cluster 2 was shown in Figure S5.

Verification of DCRS in the Guangxi Cohort and Single-Cell Dataset

To further validate the expression and predictive potential of DCRS, we collected transcriptome data from the cancer and adjacent tissues of 13 HCC patients, which we named the “Guangxi cohort” (Table S3). Tumor efficacy was evaluated based on mRECIST criteria. In the TCGA database, genes such as CLEC1B, SLC4A10, ALLC, SLC6A20, CHRDL1, and SLC34A2 showed higher expression in adjacent tissues than in tumor tissues, while genes including CRIP1, LTF, PLA2G7, CYP19A1, CHGA, NR0B1, and CYP17A1 exhibited higher expression in tumor tissues than in adjacent tissues (Figure 8A). In the Guangxi cohort, genes such as CLEC1B, SLC4A10, ALLC, CHRDL1, and SLC34A2 exhibited significantly higher expression in adjacent tissues compared to tumor tissues, while CYP19A1 showed higher expression in tumor tissues than in adjacent tissues, consistent with the findings from TCGA (Figure 8B). Additionally, in the Guangxi cohort treated with ICI and TKI combination therapy, the response group had a higher risk score than the non-response group (P=0.034) (Figure 8C), consistent with the results from GSE109211 and GSE202069. Figure 8D demonstrated the correlation of DCRS genes. We also extracted HCC tumor and adjacent samples from the single-cell dataset GSE140228 for dimension reduction, clustering, and annotation (Figure 8E–G), illustrating the expression distribution and levels of DCRS genes in various cell subgroups of the tumor microenvironment (Figure 8H). In the GSE140228 dataset, SLC4A10 showed high expression in adjacent tissues, while PLA2G7 and CRIP1 exhibited high expression in HCC, consistent with the results from TCGA and the Guangxi cohort (Figure S6). Immunohistochemical images for DCRS genes were obtained from the Human Protein Atlas database. These results were consistent with the above results (Figure S7).

Figure 8 Expression and validation of DCRS genes in Guangxi cohort. (A), Expression levels of DCRS genes in TCGA; (B), Expression levels of DCRS genes in the Guangxi cohort; (C), Boxplots of risk scores in the Guangxi cohort; (D), Correlation heatmap of DCRS genes. The “r” represents the correlation coefficient; (E-G), Dimensionality reduction, clustering, and cell annotation of the single-cell dataset GSE140228; (H), Differential expression of DCRS genes in various cells in the GSE140228 dataset.

Abbreviations: DCRS, dendritic cell-related signature; DC, dendritic cell; HSC, hematopoietic stem cell; NK, natural killer cell; TCGA, The Cancer Genome Atlas.

Notes: *P<0.05, ** P<0.01, *** P<0.001, **** P<0.0001.

Discussion

Various etiologies, the ongoing accumulation of genomic and epigenomic alterations in hepatocytes, along with alterations in the tumor microenvironment, contribute to the heterogeneity of HCC.25 Consequently, given the current rapid advancements in liver cancer treatment, there is a growing need for precision therapy tailored to HCC. DCRS has been suggested as a prognostic and therapeutic marker for colon cancer and bladder cancer.26,27 In this investigation, we identified subtypes related to DC using multiple HCC transcriptome datasets. Subsequently, we constructed and validated the DCRS, demonstrating its predictive potential for survival and its association with tumor immune infiltration. DCRS can discern HCC patients who stand to gain from chemotherapy, targeted therapy, and immunotherapy. Ultimately, consistent outcomes were derived from the HCC patients from Guangxi cohort that subjected to TKI and ICI treatments.

By integrating findings from diverse studies on risk factors, pathological features, transcriptome dysregulation, and genetic alterations, HCC can be categorized into two main groups: “Proliferation” and “Non-proliferation”.25,28–30 The “Proliferation” category is marked by poor differentiation, frequent vascular invasion, and an unfavorable prognosis, characterized by a high prevalence of TP53 mutations. Conversely, the “Non-proliferation” category is characterized by greater heterogeneity, displaying lower invasiveness, higher differentiation, and a tendency to preserve hepatocyte-like characteristics, including CTNNB1 mutations and activation of the Wnt pathway. In this study, the high-risk group, identified by DCRS, was predominantly marked by TP53 mutations (40%), aligning with the characteristics of the “Proliferation” category. Conversely, the low-risk group was primarily marked by CTNNB1 mutations (29%), resembling the characteristics of the “Non-proliferation” category. Studies suggest an antagonistic relationship between CTNNB1 mutations and TP53 mutations.31 Furthermore, integrating DCRS with clinical parameters unveiled a notably higher risk score for late-stage HCC (Stage III) compared to Stage I, with the risk score for poorly differentiated (grade 4) HCC significantly surpassing that for well-differentiated (grade 1) HCC. This implied that variations in DC infiltration characteristics may serve as a complement to the conventional molecular subtyping of HCC. As a result, differences in the abundance of TP53 and CTNNB1 mutations may form the genetic foundation for variations in clinical staging, pathological grade, immune cell infiltration, and long-term prognosis between high and low-risk groups.

HCC represents a classic inflammation-associated cancer that develops within an immunosuppressive environment, making ICIs targeting tumor immunity an attractive therapeutic approach. However, numerous clinical studies have revealed the inefficacy of PD-L1 expression levels in tumor cells for accurately reflecting tumor response.32,33 Notably, research has demonstrated the potential of the tumor immune microenvironment in predicting the response to immunotherapy. For instance, Andrew et al found that in patients with advanced HCC treated with atezolizumab in combination with bevacizumab, the T-effector signature and intratumoral CD8+ T cell density were associated with better clinical outcomes.34 The M2 subtype of tumor-associated macrophages plays a crucial role in immune suppression and tumor angiogenesis during tumor progression.35,36 Furthermore, studies have indicated that a macrophage-related signature can reflect the response to immunotherapy.37 Our investigation found that DCRS could serve as indicators of treatment responses in HCC patients undergoing TACE, TKI, and ICI therapies. Among TACE-treated patients, responders exhibited lower risk scores compared to non-responders, whereas in HCC patients subjected to TKIs or ICIs treatment, responders demonstrated higher risk scores relative to non-responders. Additionally, our results aligned with those observed in the Guangxi cohort of patients receiving TKIs and anti-PD-1 combination therapy. Immune cell infiltration analysis using xCell revealed a significant enrichment of various DC subtypes (such as pDC, iDC, and aDC) in the low-risk group. However, there were no discernible differences in the enrichment levels of immune effector cells between the high and low risk groups, including CD8+ T cells and Macrophages M1. Conversely, immunosuppressive Macrophages M2 were notably enriched in the low-risk group, potentially elucidating the suboptimal response to immunotherapy in this cohort. Studies have identified that HCC characterized by CTNNB1 mutations exhibits immune exclusion traits,38 potentially underlying the poor response to immunotherapy in the low-risk group. Beyond efficacy prediction, the study of DC enrichment has been explored as a means of advancing immunotherapy. DCs, pulsed with tumor cell lysates or transfected with a vector expressing tumor-associated antigens, have been employed to enhance the immunogenicity of secreted oncofetal antigens.39–41

In our study, the DCRS was derived from 13 genes associated with DCs. Among these, CLEC1B, CRIP1, LTF, SLC4A10, and CYP17A1 were protective factors for HCC prognosis, whereas PLA2G7, CYP19A1, CHGA, NR0B1, SLC6A20, CHRDL1, and SLC34A2 were considered risk factors. Previous research has identified CLEC1B, CRIP1, and CHGA as prognostic markers for HCC.42–44 LTF was a potential ligand targeting HCC, which could specifically bind to the asialal glycoprotein receptor. Studies used LTF as a carrier for doxorubicin to target HCC for treatment.45 Additionally, CYP17A1 and CYP19A1 were highly expressed in HCC and may serve as diagnostic markers for early-stage disease.46,47 PLA2G7 was predominantly expressed in tumor-associated macrophages and strongly correlated with poor prognosis and resistance to immunotherapy in HCC patients.48 NR0B1 enhanced HCC resistance to Sorafenib by promoting autophagy and suppressing apoptosis.49 SLC34A2 expression was upregulated in HCC cell lines, and knockdown of the SLC34A2 gene significantly inhibits HCC cell proliferation, migration, invasion, and epithelial-mesenchymal transition.50

It was also important to acknowledge the limitations of our study. One limitation was the retrospective nature of our analysis, which could introduce bias and limit the generalizability of our findings. Additionally, while we identified associations between DCRS and treatment responses, further prospective studies are needed to validate these findings and elucidate the underlying mechanisms. Moreover, our study focused primarily on transcriptomic data, and integrating other omics data, such as proteomics and metabolomics, could provide a more comprehensive understanding of the molecular landscape of HCC. Furthermore, the clinical population in this study was limited to a specific geographical region and had a relatively small sample size, which may not fully represent the global diversity of HCC patients. Future research with larger and more diverse cohorts is needed to validate the utility of DCRS as a predictive tool for HCC prognosis and treatment response.

Conclusion

This study identified genes and subtypes associated with DCs in patients with HCC. Furthermore, the DCRS model was established and validated to predict the prognosis of HCC patients, demonstrating strong predictive capability. Additionally, we assessed the responsiveness of chemotherapy, targeted therapy, and immunotherapy in DCRS risk groups. These findings have the potential to deepen our comprehension of DC characteristics and offer novel precision treatment approaches.

Abbreviations

HCC, Hepatocellular carcinoma; TKIs, Tyrosine kinase inhibitors; ICIs, Immune checkpoint inhibitors; PD-L1, Programmed death-ligand 1; TMB, Tumor mutational burden; DCs, Dendritic cells; DCRS, Dendritic cell-related signature; LASSO, Least absolute shrinkage and selection operator; TCGA, The Cancer Genome Atlas; TACE, Transarterial chemoembolization; WGCNA, Weighted gene co-expression network analysis; DEGs, Differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto; Encyclopedia of Genes and Genomes; GSEA, Gene set enrichment analysis; AIC, Akaike information criterion; OS, Overall survival.

Data Sharing Statement

The datasets used and/or analyzed during the current study are available from the corresponding author, Prof. Tao Peng, upon reasonable request.

Ethics Approval and Consent to Participate

All patient in the Guangxi cohort had signed the informed consent. The investigation had been approved by the ethics committee of Guangxi Medical University the first affiliated hospital (Approval number: 2024-E127-01). All methods in this research were carried out in accordance with Declaration of Helsinki.

Acknowledgments

The authors thank the contributors of TCGA, GSE14520, GSE76427, GSE104580, GSE109211, GSE140228 and GSE202069 database for sharing the HCC dataset on open access. We would like to thank all the HCC patients who provided specimens for this study.

Author Contributions

All authors made a significant contribution to the work reported, whether in the conception, study design, execution, acquisition of data, analysis, and interpretation, or in all these areas. They took part in drafting, revising, or critically reviewing the article, gave final approval of the version to be published, have agreed on the journal to which the article has been submitted, and agree to be accountable for all aspects of the work.

Funding

This work was supported in part by the Guangxi Research Basic Ability Improvement Project for Young and Middle-aged Teachers (2024KY0128) and the National Funded Postdoctoral Researcher Program (GZC20230583). The research was supported by the National Natural Science Foundation (82060434) and the Joint Project on Regional High-Incidence Diseases Research of Guangxi Natural Science Foundation (No.2023JJD140130). Besides, this research was partly supported by professor Tao Peng with Key Laboratory of early Prevention & Treatment for regional High Frequency Tumor (Guangxi Medical University)-Ministry of Education (grant nos. GKE2018-01, GKE2019-11, GKE-ZZ202009 and GKEZZ202109), Key R&D Plan of Qingxiu District, Nanning (NO.2020056) and Self-funded Scientific Research Project of Health Commission in Guangxi Zhuang Autonomous Region (Z20210977). The Self-raised Scientific Research Fund of the Health and Family Planning Commission of the Guangxi Zhuang Autonomous Region (grant no.Z-A20230458).

Disclosure

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA. 2021;71(3):209–249. doi:10.3322/caac.21660

2. El-Serag HB. Hepatocellular carcinoma. N Engl J Med. 2011;365(12):1118–1127. doi:10.1056/NEJMra1001683

3. Forner A, Reig M, Bruix J. Hepatocellular carcinoma. Lancet. 2018;391(10127):1301–1314. doi:10.1016/S0140-6736(18)30010-2

4. Kulik L, El-Serag HB. Epidemiology and management of hepatocellular carcinoma. Gastroenterology. 2019;156(2):477–491.e471. doi:10.1053/j.gastro.2018.08.065

5. Chen W, Zheng R, Baade PD, et al. Cancer statistics in China, 2015. CA. 2016;66(2):115–132. doi:10.3322/caac.21338

6. Xia C, Dong X, Li H, et al. Cancer statistics in China and United States, 2022: profiles, trends, and determinants. Chinese Med J. 2022;135(5):584–590. doi:10.1097/CM9.0000000000002108

7. Zihan Z, Qiulin L, Jiahua Y, et al. Analysis of the epidemiological characteristics and disease burden of malignant tumors in Guangxi, 2019. Chin J Oncol Prev Treatment. 2024;16(01):66–75.

8. Sugawara Y, Hibi T. Surgical treatment of hepatocellular carcinoma. Biosci Trends. 2021;15(3):138–141. doi:10.5582/bst.2021.01094

9. Qin S, Chan SL, Gu S, et al. Camrelizumab plus rivoceranib versus sorafenib as first-line therapy for unresectable hepatocellular carcinoma (CARES-310): a randomised, open-label, international Phase 3 study. Lancet. 2023;402(10408):1133–1146. doi:10.1016/S0140-6736(23)00961-3

10. Finn RS, Qin S, Ikeda M, Galle PR, Cheng AL. Atezolizumab plus bevacizumab in unresectable hepatocellular carcinoma. N Engl J Med. 2020;382(20):1894–1905. doi:10.1056/NEJMoa1915745

11. Rizzo A, Ricci AD, Brandi G. Trans-arterial chemoembolization plus systemic treatments for hepatocellular carcinoma: an update. J Pers Med. 2022;12(11):1788. doi:10.3390/jpm12111788

12. Sahin TK, Rizzo A, Aksoy S, Guven DC. Prognostic significance of the Royal Marsden Hospital (RMH) score in patients with cancer: a systematic review and meta-analysis. Cancers. 2024;16(10):1835. doi:10.3390/cancers16101835

13. Rizzo A, Santoni M, Mollica V, et al. Peripheral neuropathy and headache in cancer patients treated with immunotherapy and immuno-oncology combinations: the MOUSEION-02 study. Expert Opin Drug Metab Toxicol. 2021;17(12):1455–1466. doi:10.1080/17425255.2021.2029405

14. Rizzo A, Mollica V, Tateo V, et al. Hypertransaminasemia in cancer patients receiving immunotherapy and immune-based combinations: the MOUSEION-05 study. Cancer Immunol Immunother. 2023;72(6):1381–1394. doi:10.1007/s00262-023-03366-x

15. Guven DC, Sahin TK, Erul E, et al. The association between albumin levels and survival in patients treated with immune checkpoint inhibitors: a systematic review and meta-analysis. Front Mol Biosci. 2022;9:1039121. doi:10.3389/fmolb.2022.1039121

16. Greten TF, Villanueva A, Korangy F, et al. Biomarkers for immunotherapy of hepatocellular carcinoma. Nat Rev Clin Oncol. 2023;20(11):780–798. doi:10.1038/s41571-023-00816-4

17. Llovet JM, Castet F, Heikenwalder M, et al. Immunotherapies for hepatocellular carcinoma. Nat Rev Clin Oncol. 2022;19(3):151–172. doi:10.1038/s41571-021-00573-2

18. Banchereau J, Steinman R. Dendritic cells and the control of immunity. Nature. 1998;392(6673):245–252. doi:10.1038/32588

19. Takenaka MC, Quintana FJ. Tolerogenic dendritic cells. Semin Immunopathol. 2017;39(2):113–120. doi:10.1007/s00281-016-0587-8

20. Han Y, Chen Z, Yang Y, et al. Human CD14+ CTLA-4+ regulatory dendritic cells suppress T-cell response by cytotoxic T-lymphocyte antigen-4-dependent IL-10 and indoleamine-2,3-dioxygenase production in hepatocellular carcinoma. Hepatology. 2014;59(2):567–579. doi:10.1002/hep.26694

21. Obeid M, Tesniere A, Ghiringhelli F, et al. Calreticulin exposure dictates the immunogenicity of cancer cell death. Nature Med. 2007;13(1):54–61. doi:10.1038/nm1523

22. Ocadlikova D, Lecciso M, Broto JM, et al. Sunitinib exerts in vitro immunomodulatory activity on sarcomas via dendritic cells and synergizes with PD-1 blockade. Front Immunol. 2021;12:577766. doi:10.3389/fimmu.2021.577766

23. Zizzari IG, Napoletano C, Botticelli A, et al. TK inhibitor pazopanib primes DCs by downregulation of the β-Catenin pathway. Cancer Immunol Res. 2018;6(6):711–722. doi:10.1158/2326-6066.CIR-17-0594

24. Mayoux M, Roller A, Pulko V, et al. Dendritic cells dictate responses to PD-L1 blockade cancer immunotherapy. Sci, trans med. 2020;12(534). doi:10.1126/scitranslmed.aav7431

25. Rebouissou S, Nault JC. Advances in molecular classification and precision oncology in hepatocellular carcinoma. J Hepatol. 2020;72(2):215–229. doi:10.1016/j.jhep.2019.08.017

26. An B, Guo Z, Wang J, Zhang C, Zhang G, Yan L. Derivation and external validation of dendritic cell-related gene signatures for predicting prognosis and immunotherapy efficacy in bladder urothelial carcinoma. Front Immunol. 2022;13:1080947. doi:10.3389/fimmu.2022.1080947

27. Ouyang Y, Yu M, Liu T, et al. An Activated Dendritic-Cell-Related Gene Signature Indicative of Disease Prognosis and Chemotherapy and Immunotherapy Response in Colon Cancer Patients. Int J Mol Sci. 2023;24(21):15959. doi:10.3390/ijms242115959

28. Boyault S, Rickman DS, de Reyniès A, et al. Transcriptome classification of HCC is related to gene alterations and to new therapeutic targets. Hepatology. 2007;45(1):42–52. doi:10.1002/hep.21467

29. Hoshida Y, Nijman SM, Kobayashi M, et al. Integrative transcriptome analysis reveals common molecular subclasses of human hepatocellular carcinoma. Cancer Res. 2009;69(18):7385–7392. doi:10.1158/0008-5472.CAN-09-1089

30. Lee JS, Chu IS, Heo J, et al. Classification and prediction of survival in hepatocellular carcinoma by gene expression profiling. Hepatology. 2004;40(3):667–676. doi:10.1002/hep.20375

31. Schulze K, Imbeaud S, Letouzé E, et al. Exome sequencing of hepatocellular carcinomas identifies new mutational signatures and potential therapeutic targets. Nature Genet. 2015;47(5):505–511. doi:10.1038/ng.3252

32. El-Khoueiry AB, Sangro B, Yau T, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, Phase 1/2 dose escalation and expansion trial. Lancet. 2017;389(10088):2492–2502. doi:10.1016/S0140-6736(17)31046-2

33. Zhu AX, Finn RS, Edeline J, et al. Pembrolizumab in patients with advanced hepatocellular carcinoma previously treated with sorafenib (KEYNOTE-224): a non-randomised, open-label Phase 2 trial. Lancet Oncol. 2018;19(7):940–952. doi:10.1016/S1470-2045(18)30351-6

34. Zhu AX, Abbas AR, de Galarreta MR, et al. Molecular correlates of clinical response and resistance to atezolizumab in combination with bevacizumab in advanced hepatocellular carcinoma. Nature Med. 2022;28(8):1599–1611. doi:10.1038/s41591-022-01868-2

35. Biswas SK, Mantovani A. Macrophage plasticity and interaction with lymphocyte subsets: cancer as a paradigm. Nat Immunol. 2010;11(10):889–896. doi:10.1038/ni.1937

36. Wildes TJ, Dyson KA, Francis C, et al. Immune escape after adoptive T-cell therapy for malignant gliomas. Clin Cancer Res. 2020;26(21):5689–5700. doi:10.1158/1078-0432.CCR-20-1065

37. Wang T, Dai L, Shen S, et al. Comprehensive molecular analyses of a macrophage-related gene signature with regard to prognosis, immune features, and biomarkers for immunotherapy in hepatocellular carcinoma based on WGCNA and the LASSO algorithm. Front Immunol. 2022;13:843408. doi:10.3389/fimmu.2022.843408

38. Harding JJ, Nandakumar S, Armenia J, et al. Prospective genotyping of hepatocellular carcinoma: clinical implications of next-generation sequencing for matching patients to targeted and immune therapies. Clin Cancer Res. 2019;25.

39. Butterfield LH, Ribas A, Potter DM, Economou JS. Spontaneous and vaccine induced AFP-specific T cell phenotypes in subjects with AFP-positive hepatocellular cancer. Cancer Immunol Immunother. 2007;56(12):1931–1943. doi:10.1007/s00262-007-0337-9

40. Meng WS, Butterfield LH, Ribas A, et al. alpha-Fetoprotein-specific tumor immunity induced by plasmid prime-adenovirus boost genetic vaccination. Cancer Res. 2001;61(24):8782–8786.

41. Palmer DH, Midgley RS, Mirza N, et al. A Phase II study of adoptive immunotherapy using dendritic cells pulsed with tumor lysate in patients with hepatocellular carcinoma. Hepatology. 2009;49(1):124–132. doi:10.1002/hep.22626

42. Jing Q, Yuan C, Zhou C, et al. Comprehensive analysis identifies CLEC1B as a potential prognostic biomarker in hepatocellular carcinoma. Can Cell Inter. 2023;23(1):113. doi:10.1186/s12935-023-02939-1

43. Wang J, Zhou Y, Zhang D, et al. CRIP1 suppresses BBOX1-mediated carnitine metabolism to promote stemness in hepatocellular carcinoma. EMBO J. 2022;41(15):e110218. doi:10.15252/embj.2021110218

44. Wu X, Jin B, Liu X, Mao Y, Wan X, Du S. An immune-related biomarker index for predicting the effectiveness of immunotherapy and prognosis in hepatocellular carcinoma. J Cancer Res Cli