Differential analysis

Differential expression analysis was conducted using the limma package in R software. Genes with |logFC| >1 and false discovery rate (FDR) <0.01 were identified as differentially expressed genes (DEGs). We explored up- and down-regulated genes and identified differentially expressed cancer-related genes (DECRs) in HCC by intersecting them with cancer-related genes.

Screening and construction of a prognostic CR-related model

The expression levels, survival time, and survival status of the aforementioned DECRs were compiled across 370 samples. Using the “caret” R package, the samples were randomly and evenly partitioned into a training group and a testing group. Subsequently, various clinical traits of the 2 groups were examined for significant differences.

First, we performed univariate Cox regression analysis with a significant P value <0.05 to obtain candidate DECRs in the training group. Then, least absolute shrinkage and selection operator (LASSO) regression was performed on these DECRs to further remove the less relevant DECRs, and finally multivariate Cox regression was performed with the selected DECRs to build a risk score model.

The risk score for each sample is calculated as: $\text{Risk\hspace{0.17em}score}={\displaystyle {\sum }_{i=1}^{n}Coefi\times xi}$. Coef indicates the coefficient value, and x indicates the expression level of selected DECRs.

The predictive ability of the prognostic model and validation of the model

By using the “survival” R package, univariate and multivariate independent prognostic analyses were used to plot the model with other clinical traits in forest plots and examine P values to evaluate whether the model could be used as a predictor independent of other clinical traits.

Sample risk scores were ranked from low to high; risk scores, survival, and DECRs expression in the risk model were visualized, and samples were divided into high-risk and low-risk groups using the median. Receiver operating characteristic (ROC) curves at 1, 3, and 5 years were drawn to observe the predictive effect of the model on the prognosis of HCC patients, and concordance index (C-index) curves were drawn to compare the accuracy of this model in predicting the prognosis with other clinical traits.

To evaluate the prognostic value of the model in HCC patients, we analyzed the OS and progression-free survival (PFS) of the high- and low-risk samples by Kaplan-Meier analysis to observe whether there were differences between the high- and low-risk groups.

We integrated the risk formula score with demographic and clinical variables such as age, gender, stage, grade, and tumor-node-metastasis (TNM) stage to formulate a nomogram capable of predicting the 1-, 3-, and 5-year survival probabilities for HCC patients. Subsequently, we validated its accuracy by comparing the actual survival status with the predicted outcomes.

Functional enrichment analyses and gene set enrichment analysis (GSEA)

We performed Gene Ontology (GO) functional enrichment analysis on DECRs. By using the ‘org.Hs.eg.db’ R package, we converted gene names to gene id. Then, we performed enrichment analysis to obtain the gene set enrichment results. The top 6 results in BP, cellular components (CC), and molecular function (MF) were used to draw the GO circle diagram. Statistical significance was determined by P value <0.05 and FDR <0.05 during the process.

DECRs were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis. Similar to GO functional enrichment analysis, the ‘clusterProfiler’ R package was used to obtain results for gene set enrichment. The most significant 30 pathways were plotted as histograms. Statistical significance was determined by P value <0.05 and FDR <0.25 during the process.

For GSEA, we ranked DECRs in order of high to low expression, combined with channel-related gene sets downloaded from the GSEA database. The 5 most significantly enriched pathways in the high- and low-risk groups were visualized.

Tumor immune microenvironment analysis

Combined with gene expression in various immune cells and each sample, the ‘CIBERSORT’ R package was called to calculate the relative content of immune cells in these samples. Samples with insufficient accuracy were removed at P<0.05, and differential analysis was performed in the remaining samples to observe the differences in the contents of different types of immune cells between high- and low-risk groups. Then ‘GSVA’ and ‘GSEABase’ R packages were called for single sample gene set enrichment analysis (ssGSEA) to observe which immune-related functions were different in the high- and low-risk groups.

Tumor mutation burden (TMB) and tumor immune dysfunction and exclusion (TIDE)

The gene mutation data were stratified into 2 groups based on the risk score model to examine the differences in gene mutations between the high- and low-risk groups. The findings were visualized using the ‘maftools’ R package.

All samples were divided into 2 groups according to the level of TMB for survival analysis to observe whether there was a difference in OS between the 2 groups of samples, and then the samples were divided into 4 groups for survival analysis according to the combined risk score of tumor mutation load to observe whether there was an OS difference.

The expression level of the uploaded gene in each sample was used to obtain the TIDE score of each sample in http://tide.dfci.harvard.edu/, so as to compare whether there was a difference in the potential of immune escape between the high- and low-risk groups. Then, we could predict the efficacy of immunotherapy for samples at different risks.

Drug sensitivity analysis

In addition to immunotherapy, we can observe the difference in the efficacy of traditional chemotherapeutic drugs for HCC between high- and low-risk groups using the ‘oncoPredict’ R package (P<0.001) in order to guide targeted medication for patients of different risks.

Tissue sampling and clinical data

HCC and adjacent tumor tissues (>2 cm) were procured from patients who underwent radical surgical resection of HCC at the Department of Hepatobiliary Pancreatic Splenic Surgery, Affiliated Hospital of Nantong University, Jiangsu Province, China, during the period from January 2013 to June 2020. Fresh tissues were carefully aliquoted and stored at −80 °C for subsequent total protein and RNA extraction. Formalin-fixed tissues were preserved at −4 °C for the construction of tissue microarrays. A rigorous histopathological examination conducted by two experienced pathologists confirmed the identity of all specimens as HCC.

Patient-related data, encompassing age, gender, earliest diagnosed date, survival time, survival status, and TNM staging, were extracted from the hospital’s information system. Exclusion criteria were applied to eliminate patients with Class C Child-Pugh scores, positive surgical margins, preoperative interventions, history of radiotherapy, chemotherapy, targeted therapy, immunotherapy, or incomplete clinical data. Follow-up data were meticulously collected to document postoperative survival.

IHC was performed on 180 pairs of HCC and adjacent tumor tissues (see Table S1). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). All samples were acquired with informed consent from the patients, and the study protocol received approval from the Ethics Committee of the Affiliated Hospital of Nantong University (No. 2019-K021).

IHC and semiquantitative analysis

Diluted MRGBP polyclonal antibodies were evenly applied to a tissue chip containing 180 pairs of HCC tissues and adjacent tumor liver tissues. Subsequently, the chip was incubated overnight in a moisture chamber at 4 °C. The corresponding secondary antibodies were uniformly added to the tissue chip, followed by a 30-minute incubation at room temperature. We captured 3 representative images from each sample using a microscope and subjected to double-blind analysis by 2 senior pathologists, both blinded to patient clinical outcomes.

The H-score was determined based on the intensity of nuclear staining and the proportion of labeled tumor cells (21). Briefly, nuclear staining intensity was graded as 0 (no staining), 1 (weak), 2 (moderate), 3 (strong), and employed in the following formula: (% of positive cells, intensity 3 × 3) + (% of positive cells, intensity 2 × 2) + (% of positive cells, intensity 1 × 1) = H-score.

Validation of CR-based signature

The Kaplan-Meier survival curves for all groups are presented in Figure 3A-3F, revealing that patients with a high-risk score tended to exhibit lower survival probabilities and experienced earlier mortality (or metastasis) compared to those with a low-risk score. The area under the curve (AUC) values for the remaining dataset exceeded 0.65, with the exception of the 5-year AUC values, which were 0.604 for the testing group (Figure 3G-3I). These findings indicate that our risk model demonstrates a robust predictive effect on the prognosis of HCC patients.

Open in a separate window

Figure 3

Evaluation risk model. (A-C) OS of the high- and low-risk samples; (D-F) PFS of the high- and low-risk samples; (G-I) ROC curves at 1-, 3-, and 5-year; (J) comparison of ROC curves for various clinical traits to risk models; (K) comparison of C-index curves for various clinical traits to risk models; (L) univariate independent prognostic analysis to assess the risk model; (M) multivariate independent prognostic analysis to assess the risk model; (N) nomogram containing various clinical traits and risk scores. ***, P<0.001; (O) calibration images of predicted versus actual values. OS, overall survival; PFS, progression-free survival; ROC, receiver operating characteristic; C-index, concordance index; AUC, area under the curve; CI, confidence interval.

In the 3-year ROC curves, the AUC values for risk scores surpassed those for other clinical traits, suggesting that the utilization of risk scores provides a more accurate prediction of the survival of HCC patients compared to other clinical traits (Figure 3J,3K).

Both univariate and multivariate independent prognostic analyses showed P<0.05, indicating that the model can be used as a predictive tool independent of other clinical traits (Figure 3L,3M).

To predict the 1-, 3-, and 5-year survival probabilities for each sample, we constructed a nomogram containing various clinical traits and risk scores (Figure 3N). A higher score means a lower probability of survival. For example, a patient with HCC who scored 392 had a 1-year survival rate of 0.929, a 3-year survival rate of 0.861, and a 5-year survival rate of 0.81. Calibration images showed that nomograms predicted 1-, 3-, and 5-year OS were in good agreement with actual values and accuracy (Figure 3O).

In summary, this risk model has been validated as a tool to predict prognosis.

Functional enrichment analyses and GSEA

The top 6 with the highest significance among the GO enrichment analysis results were selected and integrated into a circle diagram (Figure 4A). As can be seen in the figure 4A, in BP, these DECRs were mainly involved in these biological functions such as nuclear division, organelle fission, mitotic nuclear division, chromosome segregation, nuclear chromosome segregation, and sister chromatid segregation; in CC, cell components such as chromosomal region, chromosome, centromeric region, condensed chromosome, the CDC45/RecJ, MCM, GINS (CMG) complex, and DNA replication preinitiation complex were involved; in MF, these risk DEGs were mainly involved in MFs such as single-stranded DNA helicase activity, tubulin binding, microtubule binding, DNA helicase activity, ATP-dependent activity acting on DNA, and catalytic activity acting on DNA.

Open in a separate window

Figure 4

Functional enrichment analysis of differential genes. (A) GO circle diagram of DEGs; (B) histogram of KEGG functional enrichment analysis; (C) enriched gene sets in the high-risk group; (D) enriched gene sets in the low-risk group. GO, Gene Ontology; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

KEGG pathway analysis showed these DECRs were significantly associated with cell cycle pathways, DNA replication, extracellular matrix (ECM)-receptor interaction, and protein digestion and absorption (Figure 4B).

Via GSEA, we noted enrichments in cell cycle, cytokine receptor interaction, ECM receptor interaction, hematopoietic cell lineage, and neuroactive ligand receptor interaction in the high-risk group. Conversely, the low-risk group exhibited enrichments in drug metabolism cytochrome P450, fatty acid metabolism, glycine serine and threonine metabolism, primary bile acid biosynthesis, and retinol metabolism (Figure 4C,4D).

Tumor immune microenvironment analysis

After selecting eligible samples, their relative immune cell content was visualized and analyzed for differences. It could be observed that monocytes were significantly more frequent in low-risk samples than in high-risk samples, whereas macrophages M0 cells were significantly more frequent in high-risk samples than in low-risk samples (Figure 5A). Among various immune-related functions, B cells, cytolytic activity, mast cells, neutrophils, NK cells, T helper cells, type I interferon (IFN) response, and type II IFN response were significantly active in the low-risk group compared with the high-risk group; activated dendritic cells (aDCs), antigen-presenting cell (APC) co stimulation, macrophages, major histocompatibility complex (MHC) class I, T follicular helper (Tfh), and tumor-infiltrating lymphocyte (TIL) was significantly active in the high-risk group compared with the low-risk group (Figure 5B).

Open in a separate window

Figure 5

Tumor immune microenvironment analysis. (A) Relative immune cell content of eligible samples; (B) comparison of immune-related functions in high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001. aDCs, activated dendritic cells; APC, antigen-presenting cell; CCR, C-C chemokine receptor; DCs, dendritic cells; HLA, human leukocyte antigen; iDCs, immature dendritic cells; MHC, major histocompatibility complex; pDCs, plasmacytoid dendritic cells; Tfh, T follicular helper; TIL, tumor-infiltrating lymphocyte; IFN, interferon; NK, natural killer.

TMB and TIDE

TMB analysis of samples from the high- and low-risk groups showed that the gene tp53 had the most mutations in samples from the high-risk group whereas the gene CTNNB1 had the most mutations in samples from the low-risk group. Missense mutations were the most frequent mutation type in most of the mutated genes in both groups (Figure 6A).

Open in a separate window

Figure 6

TMB and TIDE. (A) TMB analysis; (B) survival probability of samples in high- and low-TMB groups; (C) survival probability of samples in high TMB + high risk group, high TMB + low risk group, low TMB + high risk group and low TMB + low risk group; (D) TIDE scores for high- and low-risk groups. TMB, tumor mutation burden; TIDE, tumor immune dysfunction and exclusion. ***, P<0.001.

At the same time point, significantly fewer samples survived in the high TMB group than in the low TMB group. In the 4 groups that combined the TMB with the risk model score to divide the samples, it could be observed that the low TMB + low risk group was located at the top of the graph whereas the high TMB + high risk group was located at the bottom of the graph meaning that they had the best and worst prognosis, respectively (Figure 6B,6C).

A notable distinction in TIDE scores was observed between the high- and low-risk groups, with the low-risk group exhibiting lower TIDE scores. This suggests that individuals with low-risk HCC may experience enhanced effectiveness with immunotherapy (Figure 6D).

Drug sensitivity analysis

The results showed that PF-4708671, JQ1, ribociclib, and SB505124 could obtain better efficacy in the treatment of high-risk group samples, whereas taselisib, ipatasertib, GDC0810, cediranib, and 5-fluorouracil could obtain better efficacy in the treatment of low-risk group samples. These findings are conducive to the development of accurate treatment plans for high- and low-risk groups, respectively (Figure 7).

Open in a separate window

Figure 7

Prediction of differential chemotherapy drug sensitivity between high- and low-risk groups.

Funding: This study was supported by the National Natural Science Foundation of China (Nos. 81871927, 81070360), 2021 Changchun University of Chinese Medicine School-level Clinical Practice Teaching Reform Special Research Project (No. XJLCSJ202146), Scientific Research Project of Nantong Health Commission (No. MS12020024), Jiangsu Province Graduate Research Innovation Plan Project (No. KYCX23_3423), and the Nantong Hepatobiliary and Pancreatic Surgery Disease Research Center Construction Project (No. HS2015001).