We obtained HCC patients’ RNA sequencing (RNA-seq) expression data, genomic mutation data, and accompanying clinical data from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/repository). Similarly, the RNA-seq expression data and matched clinicopathologic information were also downloaded from the GSE14520 database (16) of the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/gds/) and from the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/projects/LIRI-JP). We utilized the transcriptome data and clinical data of the IMvigor210 dataset of patients with metastatic urothelial cancer treated with anti-PD-L1 drugs obtained from the IMvigor210CoreBiologies R package to verify whether the risk model we established could predict the effectiveness of immunotherapy. The flowchart of the present study design is shown in Figure 1.

The TIMER algorithm was utilized to estimate the relative proportion of infiltrating immune cells in HCC samples in each cohort (17). We used WGCNA which can convert coexpression correlation into connection weights or topological overlap values (18), to identify coexpressed genes in macrophages. Network type was set as the “unsigned” type. We used standard deviation (SD > 50%) to screen highly variable genes in the WGCNA expression data. The topological overlap matrix was employed to determine the connectivity and dissimilarity of the coexpression network established with the appropriate soft thresholding parameter (Figure 1). The dynamic hybrid cut method (a bottom-up algorithm) was used to identify coexpressed gene modules. A hierarchical clustering tree was established using dynamic hybrid cutting. Each leaf on the tree represents a single gene, and genes with similar expression data or similar functions are close together and form a branch of the tree, representing a gene module. We used Pearson’s test to calculate the correlation between module eigengenes (MEs) and macrophages. When P <0.05, the module was considered to be significantly related to macrophages. In this way, a set of genes related to the proportion of macrophages with similar functions was identified in the same module. Subsequently, we took the intersection of the genes of the modules most related to macrophages in the three cohort datasets and displayed them in the form of a Venn diagram.

Based on the expression of macrophage-related intersection genes, unsupervised hierarchical clustering was applied to classify HCC patients into the optimal number of clusters using the “ConsensusClusterPlus” package of R (19).

Tumor mutational burden (TMB) in each tumor sample refers to the number of mutated bases per million bases, which include missense mutations, nonsense mutations, frameshift mutations, and so on. We computed the TMB values from the number of variants out of the total length of the human exons (38 million) in each sample by using Perl scripts. Waterfall plots were generated using the “maftools” R package to assess the number of somatic point mutations in each HCC sample and to illustrate the relationship between TMBs and risk groups (20). We calculated the copy number variation (CNV) frequencies and displayed the above results in a lollipop chart. The “RCircos” package of R software was utilized to visualize the location of these genes on the chromosomes.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used to assess the biological roles of the prognostic candidates by the “clusterProfiler” R package. Gene set enrichment analysis (GSEA) software (version 4.1.0) was used to compare the biological processes that were significantly different between the low- and high-risk groups.

In order to investigate the difference in the TME among macrophage-related clusters, the stromal score, immune score, estimate score, and tumor purity were estimated using the “ESTIMATE” R package based on the results of single-sample gene set enrichment analysis (ssGSEA) (21). ssGSEA with the “gsva” package was implemented to calculate the activity scores of the activities of immune-related pathways (22). Immune infiltration and functions were compared between the different groups by using the two-sample Wilcoxon test. We further compared the expression of common immune checkpoint inhibitor (ICI) and human leukocyte antigen (HLA) genes between different clusters. Additionally, we acquired an overview of the immune subtypes of HCC patients from the TGGA database from UCSC-Xena (https://xenabrowser.net/datapages/) (23). We made comparisons of the immune subtypes between different risk groups based on macrophage-related clusters using the package “RColorBrewer.”

The immunophenoscore (IPS) of HCC samples was obtained from the LIHC project of The Cancer Immunome Atlas (TCIA, https://tcia.at/), which can predict the response to immunotherapies including CTLA4 and PD-1 blockers (24). The IPS score was normalized to a range of 0 to 10, where a higher IPS score represents higher immune reactivity. The tumor immune dysfunction and exclusion (TIDE) was used to predict the potential immune checkpoint blockade responses in HCC. A lower TIDE score represents a better response to immunotherapy. In addition, we utilized the IMvigor210 dataset to validate the links between the risk signature and immunotherapy.

The TISMO (tismo.cistrome.org) database was used to compare gene expression levels across groups of different responses to immune checkpoint blockers (ICBs) in syngeneic mouse models.

Based on the RNA-seq results, differentially expressed genes (DEGs) between cluster 1 and cluster 2 were identified. Then, we used univariate Cox regression analysis to identify the genes with good predictive ability for prognosis. To further narrow down the candidate prognosis-related genes and synthetically estimate the significance values, the LASSO Cox regression algorithm (R package “glmnet”) was utilized to identify the variation in regression coefficients of the prognostic genes and select the optimal and minimal criteria of the penalization parameter (λ) using 10-fold cross-validation. Immunohistochemistry data from the Human Protein Atlas (HPA) database confirmed the expression of these model genes in HCC. All HCC patients were grouped into high- and low-risk groups by the median risk score value. The difference in expression between the identified genes by the LASSO Cox regression algorithm and the distribution pattern between clinicopathological characteristics and risk groups are displayed in the form of heatmaps using the R package “pheatmap.” Principal component analysis (PCA) and t-distributed random neighborhood embedding (t-SNE) analysis with the R package “Rtsne” were applied to gauge the discriminative ability of the predictive model. The area under the time-dependent receiver-operating characteristic curve (AUROC) was used to appraise the predictive ability of the risk groups identified above. A risk curve was drawn to explore the difference in survival status between different risk groups of patients. Then, the accuracy of the risk score model was validated in the GSE14520 and ICGC datasets using the same method.

To evaluate whether risk score is an independent risk factor affecting the survival of HCC patients, univariate and multivariate analyses were conducted on the TCGA, GSE14520, and ICGC datasets. Similarly, decision curve analysis (DCA) was carried out to determine the clinical application value of the risk score model by calculating the net benefits at each risk threshold (25). A nomogram integrating sex, grade, age, stage, and risk signature for survival prediction was then established using the “regplot” package. The calibration curve was drawn to evaluate the predictive accuracy of the nomogram.

In view of the lack of biomarkers that can accurately predict chemotherapeutic drug sensitivity in HCC patients, we conducted drug susceptibility analyses through the “pRRophetic” and “ggplot2” packages to compare the half-maximal inhibitory concentration (IC₅₀) of various chemotherapeutic drugs for HCC between the high- and low-risk groups using the Wilcoxon signed-rank test (26).

All statistical analyses in our study were executed by R software (version 4.1.0, https://www.r-project.org/). For comparison of continuous variables between two groups, the independent Student’s t-test was performed, and the Wilcoxon rank-sum test was utilized to compare non-normally distributed variables. The chi-squared test was used for the comparison of categorical variable data between two groups. The P-value for statistically significant differences was set to 0.05, unless otherwise stated in the manuscript.

In the TCGA-LIHC cohort dataset, among the 10 modules, the yellow module was highly correlated with macrophages (R² = 0.51, P = 1e−26). In the GSE14520 cohort dataset, among the eight modules, the yellow module was highly correlated with macrophages (R² = 0.51, P = 3e−26). In the ICGC cohort dataset, among the nine modules, the turquoise module was highly correlated with macrophages (R² = 0.58, P = 1e−22). Given that our study focused on macrophages, the yellow module from the TCGA-LIHC dataset and the GSE14520 dataset and the turquoise module from the ICGC dataset were identified as the hub module, as shown in Figure 2A.

Subsequently, 174 macrophage coexpressed genes were identified from the intersection of the above modules (Figure 2). We compared the expression levels of all 174 macrophage coexpressed genes between tumor tissues and normal samples, and we identified 69 differentially expressed macrophage coexpressed genes (|log₂ FC|  ≥ 0.585, all P < 0.05), as shown in Supplementary Figure S1. Then, the potential interactions among these differentially expressed macrophage coexpressed genes were analyzed by the PPI network (Supplementary Figure S1). The correlation network containing all differentially expressed macrophage coexpressed genes is presented in Supplementary Figure S1. GO enrichment analysis and KEGG pathway analysis were then applied to the genes lying in the intersection of the above three cohort datasets (Figure 2K).

To further explore the heterogeneity of macrophages between tumors in HCC, based on the expression level of 174 macrophage coexpresssed genes, unsupervised consensus analysis was employed to divide these 174 macrophage coexpression genes in TCGA-LIHC into two different subtypes with different molecular and clinical characteristics, including 226 samples in macrophage-related cluster 1 and 144 samples in macrophage-related cluster 2 (Figures 3A–C). 3D PCA showed that the two clusters could be well-distinguished by the macrophage-related coexpressed genes (Figure 3D). The heatmap indicated the distribution of clinicopathological characteristics in cluster 1 and cluster 2. We found significant differences in tumor stage status and pathological grade in the two clusters (Figure 3E). Furthermore, we found that patients in cluster 1 had significantly better OS than those in cluster 2 (Figure 3F). GSVA enrichment analysis showed that there were differences in the enriched pathways between cluster 1 and cluster 2 (Figure 3G). In total, 2,216 DEGs between the different clusters were identified. GO enrichment and KEGG pathway analyses were then run on the above DEGs, as shown in Figures 3H–K.

To verify the reliability of the above cluster scheme and compare the tumor immune microenvironment and activities of immune-related pathways between the two groups, we calculated the immune score, estimate score, stromal score, and tumor purity based on the gene expression profile of each HCC sample using the ESTIMATE algorithm. The immune score, stromal score, and estimate score in cluster 1 were significantly lower than those in cluster 2 (Figure 4), while tumor purity in cluster 1 was significantly higher than that in cluster 2 (Figure 4). These results suggested that the above two clusters had completely different TME infiltration characteristics. There were marked differences in immune cell-related functions between the two clusters. The type II IFN response was enriched in cluster 1, but others were enriched in cluster 2. Finally, we compared the expression of ICI genes and HLA genes between clusters. Immune checkpoint genes (HAVCR2, VSIR, PDCD1, CTLA4, CD276, CD274, LAG3, NRP1, TIGHT, BTNL2, IDO1, TNFRSF14, and VTCN1) were highly expressed in cluster 2, while the HLA-related genes were highly expressed in cluster 2 (Figure 4D). To further test the ability of the macrophage-associated clusters to predict the ICI response, IPS analysis was carried out to determine the immunotherapeutic sensitivity of HCC patients, respectively. As shown in Figure 4F, the IPS–PD1 blocker score, the IPS–CTLA4 score, and the PD1 blocker score in cluster 2 were significantly higher than those in cluster 1, which indicated a more immunogenic phenotype in cluster 2, so HCC patients in cluster 2 might benefit from immunotherapy.

TMB, also known as non-synonymous variation, is strongly associated with immune cell infiltration and immune response (27). Somatic mutation data of HCC patients were obtained and TMB scores were computed. The waterfall plot showed mutations of the macrophage-related DEGs (Figure 4). We divided the HCC patients into the low-TMB and high-TMB groups around the optimal cutoff value. The Kaplan–Meier curve indicated that the prognosis of patients with high TMB was poorer than that of patients with low TMB. In addition, given the presence of a potential correlation between macrophage-related clusters and TMB, we performed the stratification analysis and found that the cluster combined with the TMB risk group can better predict the prognosis of HCC patients (Figure 4K).

To better apply these subtypes to the clinical treatment of HCC and calculate the specific risk score of each HCC patient, subsequently, we compared the DEGs between the two clusters and established a specific risk scoring model based on macrophage-related clusters (Figure 5). Subsequently, based on the obtained gene expression profiles, an 18-gene risk score model based on macrophage-related clusters was established for each patient based on the personalized gene level using LASSO Cox regression analyses (Supplementary Figure S2), named the macrophage-related risk score signature (MRS) (Figures 5B, C). We utilized the Kaplan–Meier curve analysis to compare the effects of high and low expression of the above-selected genes on long-term survival (Supplementary Figure S3). Kaplan–Meier analysis revealed that patients in the high-risk group had significantly poorer OS than those in the low-risk group, as shown in Figures 5D–G. PCA and t-SNE analysis indicated that the MRS model showed good discrimination (Figures 5H–M). Risk curve analysis indicated that the high-risk group had higher mortality and shorter survival time (on the right side of the dotted line). The number of patients in the high-risk group grew and death events increased with increasing risk score (Figures 5N–P). The AUCs for 1-, 3-, and 5-year OS were 0.808, 0.749, and 0.746, respectively, in the TCGA cohort (Figure 5Q). Consistent with the results in the TCGA cohort, the AUC values for 1-, 3-, and 5-year OS were 0.661, 0.681, and 0.667, in the GSE14520 cohort, respectively (Figure 5R). The AUC values for 1-, 3-, and 5-year OS were 0.791, 0.780, and 0.591 in the ICGC cohort, respectively (Figure 5S). Furthermore, the AUC value for MRS based on macrophage-related clusters was significantly higher than those for age, sex, tumor stage, and pathological stage in the TCGA-LIHC cohort. Meanwhile, the MRS model revealed a good predictive ability that was not inferior to the predictive ability of other clinicopathological factors in the GSE14520 and ICGC cohorts (Figures 5T–V).

Immunohistochemical results from the HPA database were utilized to further assess the expression of genes from the MRS model in HCC. Compared with adjacent normal liver tissues, the protein expression levels of genes (PON1, KLRB1, ANXA10) in HCC tissues decreased significantly, while the protein expression levels of IL15RA, HILPDA, HAVCR1, G6PD, CDCA8, and CBX2 in HCC tissues increased. In addition, these data showed significant differences in PON1, KLRB1, ANXA10, IL15RA, HILPDA, HAVCR1, G6PD, CDCA8, and CBX2 between HCC and normal liver tissues, and this trend was the same as in the model (Supplementary Figure S4).

We found that there were significant correlations between risk score and pathologic grade, T status, and tumor stage in the TCGA-LIHC cohort (Figure 6, P < 0.05). There were positive correlations between risk score and AFP level, cirrhosis, tumor size, and tumor stage in the GSE14520 cohort (Supplementary Figures S5A–G). The risk score in stage IV was significantly higher than that in stages I and II, and the risk score in the death group was significantly higher than that in the alive group in the ICGC cohort (Supplementary Figures S5H–J). In addition, univariate and multivariate Cox regression analyses indicated that the MRS was an independent factor predicting survival, as shown in Figures 7A–F. The decision curve analysis indicated that the MRS had a higher clinical net benefit than other clinicopathological characteristics in the TCGA-LIHC cohort (Figure 7G). Subsequently, based on the results of the stepwise Cox regression model, we further constructed a clinically adaptable nomogram with the MRS and other clinicopathological characteristics to provide a visual way to predict the 1-, 3-, and 5-year survival with HCC (Figure 7H). Our nomogram exhibited better accuracy in predicting both short- and long-term survival. The calibration plot of the nomogram showed excellent concordance between the prediction by the nomogram and the actual observation probabilities (Figures 7I–K). The above findings suggest that the nomogram we established has good prognostic value for patients with HCC.

To further analyze the differences in the gene functions and involved pathways between the subgroups classified by MRS, 220 DEGs between the low- and high-risk groups were identified. GO enrichment analysis and KEGG pathway analysis were then performed on these DEGs. GO pathway enrichment analysis revealed that the DEGs were mainly concentrated in “mitotic nuclear division,” “mitotic sister chromatid segregation,” “sister chromatid segregation,” and “nuclear division” in the biological process category. The KEGG analysis results showed that the DEGs were mainly enriched in pathways associated with “complement and coagulation cascades,” “cell cycle,” “retinol metabolism,” “drug metabolism-cytochrome P450”, and “metabolism of xenobiotics by cytochrome P450” (Figure 8).

GSEA was further performed to complement and validate the functional annotation of KEGG and GO. KEGG enrichment analysis indicated that the most enriched pathways in the high-risk group were “cell cycle,” “DNA replication,” “ECM receptor interaction,” and “neuroactive ligand-receptor interaction.” In contrast, “fatty acid metabolism,” “retinol metabolism,” and “drug metabolism cytochrome P450” were enriched in the low-risk groups. Additionally, GO enrichment analysis indicated that the most enriched biological processes in the high-risk group were closely linked to “humoral immune response mediated by circulating immunoglobulin,” “phagocytosis recognition,” “B-cell receptor signaling pathway,” and “positive regulation of B-cell activation” (Figure 9).

A CNV is a DNA fragment with copy sizes ranging from 1 kB to 1 MB in the human genome that is related to activation of oncogenes or inactivation of tumor suppressor genes and heterogeneity of the genome and molecular phenotype, further leading to the tumorigenesis and progression of tumors (28, 29). We analyzed the frequencies of genetic amplification and deletion of selected genes from the MRS. The above results revealed that S100A9, CBX2, STC2, G6PD, and PFN2 had a higher frequency of gain-of-function mutations in HCC, while CDCA8, NEIL3, IL15RA, SPP1, and KLRB1 had a higher frequency of loss-of-function mutations (Figures 10A, B). In view of the important role that TMB may play in clinical practice, we attempted to understand the intrinsic association between TMB and MRS. The results revealed that TMB was higher in the high-risk group and that the risk score was positively correlated with the TMB score in HCC (r=0.13, P=0.018, Figures 10C, D). Furthermore, it appeared that there was a combined influence of TMB and MRS on survival outcomes in patients with HCC (Figure 10E). Waterfall plots demonstrated the mutation differences of the top 20 genes between different risk groups in HCC (Figures 10F, G). Consistent with the TMB score, patients in the high-risk group had a higher mutation frequency. The mutation frequency of TP53 in the high-risk group was significantly higher than that in the low-risk group, while the mutation frequency of AXIN1 in the low-risk group was significantly higher, and it was mainly composed of frameshift deletion and nonsense mutation.

We compared the proportions of immune subtypes of HCC between different risk groups. The results indicated that there were significant differences in immunophenotyping between different groups (Figures 11A, B). Immune-related functions, including cytolytic activity, type I IFN response, and type II IFN response, were enriched in the low-risk group, and MHC class I responses were enriched in the high-risk group. We found that the expression of ICI genes in the high-risk group was generally higher than that in the low-risk group, except for TDO2 (Figures 11C, D). Additionally, the results revealed the correlation between MRS and ICI gene expression, as shown in Supplementary Figure S6.

We further evaluated whether the MRS could serve as an immunotherapy predictor for HCC patients. Based on the TIDE algorithm, we observed that there was a negative correlation between TIDE score and risk score, and the TIDE score in the high-risk groups was significantly lower than that in the low-risk group (Supplementary Figure S7). The percentage of responses in the high-risk group was higher than that in the low-risk group (Figures 11E–G). In the IMvigor210 cohort, 348 patients with metastatic urothelial cancer treated with anti-PD-L1 drugs were segmented into high- and low-risk groups, although no significant difference in OS was detected (P = 0.43), and the high-risk group showed a better survival trend (Figure 11K). There was no significant difference in the proportion of inflammatory immune subtypes or tumor-infiltrating immune cells expressing PD-L1 between the high- and low-risk groups (Figures 11H, I). However, compared with the low-risk groups, we found that the tumor tissue samples of the high-risk group had a higher proportion of tumor cells expressing PD-L1 (Figure 11J). In addition, we found that neoantigen count was positively correlated with risk score (Supplementary Figure S7), suggesting that a higher level of immune events might exist in patients with higher risk scores (R = 0.31, P = 6e−17). Compared with the low-risk group, the high-risk group had a higher percentage of complete response/partial response (Figures 11L, M). Analysis of two cell line data (BNL-MEA) from the TISMO database showed upregulated NEIL3, CBX2, and CDCA8 related to a poor response for ICB and downregulated NEIL3, CBX2, and CDCA8 related to a better response for ICB when compared with maintaining baseline expression level. In addition, the data revealed the downregulated TTK, S100A9, PON1, PFN2, and GPC1 related to a poor response for ICB and upregulated TTK, S100A9, PON1, PFN2, and GPC1 related to a better response for ICB when compared with maintaining baseline expression level (Supplementary Figure S8).

To study the possible application of MRS in the personalized treatment of HCC, we assessed the IC₅₀ values of several chemotherapy agents between the different MRS signature groups. The results of drug susceptibility indicated that the high-risk group had lower IC₅₀ values of gemcitabine, bleomycin, cisplatin, doxorubicin, mitomycin C, and paclitaxel than the low-risk group, which implied that patients with high risk may benefit more from the above chemotherapies, while the IC₅₀ values of sorafenib, rapamycin, bosutinib, dasatinib, docetaxel, and metformin were lower in the low-risk group, which suggested that patients with low risk may benefit more from the above chemotherapies (Figure 12).