In this study, a total of 365 HCC patients, including complete RNA transcriptome data (FPKM normalized) and clinical characteristics (including age, gender, tumor grade, T stage, N stage, and M stage), were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov) as TCGA cohort. HCC patients from a local HCC cohort (n = 115), including clinical characteristics and antigen-presenting related gene (ARG) expressions, were enrolled at the First Affiliated Hospital of Wenzhou Medical University (FAHWMU) (Wenzhou, China). The collection of this cohort was reviewed and approved by the human research ethics committee of the FAHWMU. All patients/participants provided their written informed consent to participate in this study. The 17 ARGs (PNRC1, CTSS, CTSF, HLA-DRA, HLA-DRB1, ATF3, RFX3, USP11, USP34, PSMA3, PSMA4, IRF1, CD207, CLEC4A, CLEC5A, CD1A, and CD1C) used in this study were gained from MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb).

Among the clinical characteristics in the FAHWMU cohort, the laboratory variables were taken from the results of a test closest to the date of surgery, including hepatitis B, alpha-fetoprotein (AFP), carcinoembryonic antigen (CEA), and carbohydrate antigen 19-9 (CA199). The histopathological variables (including tumor size, lymph node invasion, vascular invasion, and perineural invasion) were also included. Based on the 8th edition of the American Joint Committee on Cancer (AJCC) Staging Manual, the TNM stage for each HCC patient was obtained. In addition, age and gender were included as the demographic characteristics.

The total RNA from the liver tissues of the FAHWMU cohort was extracted using TRIzol reagent. The mRNA was then reverse transcribed into cDNA using ribo SCRIPTTM reverse transcription kit. The expression levels of mRNA were calibrated with glyceraldehyde-3-phosphate dehydrogenase (GAPDH). SYBR Green master mix was added, and real-time PCR was carried out using a 7500 rapid quantitative PCR system (Applied Biosystems, Foster City, CA, USA). The CT value of each well was recorded, and the relative quantification of the amplified products was performed using the 2^−ΔCt method.

To initially analyze the related functionality of ARGs, the Gene Ontology (GO) enrichment analysis was performed according to the Enrich database (http://amp.pharm.mssm.edu/Enrichr/) (21). In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, based on the DAVID database (https://david.ncifcrf.gov/), was utilized to reveal the ARG-related potential signaling pathways. Through the Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org) dataset, the response of HCC patients to potential chemotherapeutic drugs was predicted. The sensitivity of chemotherapeutic drugs was calculated by the half-maximal inhibitory concentration (IC50) method, and IC50 was estimated using the R package “pRRophetic” (22). In targeting correlation analyses with respect to immune checkpoints, the Wilcoxon test and Pearson’s correlation analysis were applied to calculate the relationship between different groups.

Raw data on tumor mutational burden (TMB) for HCC patients were enrolled from TCGA database (https://portal.gdc.cancer.gov/). With the aid of the “mattool” R package, upstream analysis of whole-genome sequencing and whole-exome sequencing data was performed (23). For systematic analysis of TMB indicators in individual HCC patients, somatic mutation analysis was used, with the “getsamplesummary” function and the “getgenesumprice” function retrieving sample information and gene information, respectively (24).

The consensus clustering method was used to identify the APC-related subtypes in this study (25). In the consensus clustering, the resampling method was applied to reduce the error. Meanwhile, the number of cluster k values was specified to calculate the rationality under different cluster numbers (26). The k value with the smallest cumulative distribution function (CDF) slope was considered the optimal value for separating clusters when the consensus index value ranged from 0.1 to 0.9 (27).

The relative ESTIMATE score, Immune score, and Stromal score for individual HCC patients were estimated using the Estimation of STromal and Immune cells in MAlignant Tumor tissues using the Expression data (ESTIMATE) algorithm (28). Meanwhile, the proportion of 22 tumor-infiltrating immune cells (TICs) was calculated through the CIBERSOFT algorithm, which is a gene-based deconvolution algorithm that infers 22 human TIC types and quantifies the relative fraction of each cell type (29). The details for genes, which represent 22 TICs, are presented in Table S1. To assess the activity of the programmed death-ligand 1 (PD-L1) antibody in patients, the IMvigor210 trial was also applied (30).

Kaplan–Meier (K-M) method was applied to construct K-M survival curves (31). Meanwhile, the log-rank test was used to validate the accuracy of the K-M survival curve (32). The receiver operating characteristic (ROC) curve, based on the R package “timeROC,” was generated to evaluate the correlation between indices and overall survival (OS) (33). The univariate Cox regression analysis was employed to further evaluate the clinical prognostic value. The predictors including clinical characteristics and risk scores were integrated by the nomogram. Through the Cox regression method, the calibration curve at the 1st, 2nd, and 3rd years and the PH hypothesis analysis were constructed to verify the accuracy of the nomogram Table 1 (34).

The overview of the study design is shown in Figure 1. As presented in Figure 2A, the GO function enrichment analysis revealed the ARG-enriched significant cellular functions. Except for the functions related to antigen presentation, there were also a large number of cellular functions closely related to various pathways of the immune response as well as stress response, protein ubiquitination, and interferon signature pathways. The protein-to-protein interaction (PPI) relationships of ARGs were displayed in the PPI network (Figure 2B). It was found that CD1C, IRF1, and CLEC4A had the most protein interactions with other ARGs. As indicated in Figure 2C, the ARG expressions were examined between adjacent non-tumorous samples and HCC samples. Our results showed that PNRC1, ATF3, CD207, and CD1C were obviously reduced in the HCC samples compared to non-tumorous samples, whereas RFX3, USP11, USP34, PSMA3, PSMA4, CLEC4A, and CD1A were increased in HCC samples. The correlation network of ARGs was presented in Figure 2D.

Based on the CDF method with an optimal k value of 2, HCC patients in TCGA cohort were classified into 2 APC-related subtypes via consensus clustering analysis (Figures 3A, B, Figure S1). As shown in Figure 3C, the K-M survival curves showed that patients in the APC cluster A had a significantly lower OS compared to APC cluster B (Figure 3C, p = 0.017). The complex heatmap displayed the distribution of clinical characteristics and ARG expressions between APC clusters (Figure 3D). Most ARGs showed higher expression levels in APC cluster A in comparison with APC cluster B. As shown in Figure 4A, the relative scale of a fraction of plasma cells, T-cell CD4 memory activated, T-cell follicular helper, and T-cell regulatory (Tregs) was significantly downregulated in APC cluster B compared to APC cluster A (p < 0.01). By contrast, monocytes, mast cell resting, and mast cells activated were remarkably upregulated in APC cluster B (p < 0.01). Moreover, the TME scores (Immune score, Stromal score, and ESTIMATE score) were obviously lower in APC cluster B (p < 0.01). For TME scores, lower Immune score and Stromal score indicated a lower relative content of stromal cells and immunocytes in TME, while a lower ESTIMATE score indicated the lower aggregation of Stromal and Immune scores in TME. The associations between APC clusters and immune checkpoints were next analyzed. All immune checkpoints including CTLA4, FOXP3, IL10, and PD-L1 were obviously downregulated in APC cluster B (Figure 4B, p < 0.05), implying that the immune system of patients with APC cluster A may be suppressed. In addition, a lower TMB value was shown in APC cluster B (Figure 4C, p = 0.036), indicating that patients in APC cluster A may benefit from immunotherapy.

With the use of the Wilcoxon test, a total of 980 significantly differentially expressed genes (DEGs) were screened in APC clusters as the potential genes targeting APCs. Subsequently, the consensus clustering was performed again based on the genome matrix for 980 DEGs. Based on the CDF method, the optimal k value was selected to divide HCC patients into 3 main APC-related genomic subtypes (Figures 5A, B, Figure S1). The K-M survival curve was used to analyze the OS among genomic clusters. Our data indicated that genomic cluster C had lower OS compared to genomic clusters A and B (Figure 5C, p = 0.026). We additionally analyzed the scale of fraction for 22 main TICs and ESTIMATE scores among genomic subtypes (Figure 5D). The relations between clinical characteristics and DEGs expressions in different genomic subtypes are shown in the heatmap (Figure 5E). The expression of most DEGs tended to increase as APC-related clustering from A to C. The GO enrichment analysis showed the functions of immune responses, multiple metabolic responses, and various basic physiological activities differed significantly among genomic subtypes (Figures 6A, B, p < 0.05).

Based on the Lasso-Cox regression analysis of DEGs, a total of 12 hub prognostic-related genes were screened to generate the prognostic model. The risk score under the prognostic model was next calculated according to the gene expressions and optimal coefficients.

Risk score = EIF2S1 * 0.0087

+YBX1 * 0.0011

+OLA1 * 0.0022

+TRMT10C * 0.0032

+KIF20A * 0.0075

+SF3B4 * 0.0047

+ETF1 * 0.0024

+CDCA8 * 0.013

+YARS1 * 0.029

+DLAT * 0.0034

+TIMM23 * 0.0056

+GDI2 * 0.0020

All HCC patients in TCGA cohort were divided into the high-risk group (n = 149) and low-risk group (n = 216) according to the optimal risk score. The expression levels of 12 genes between risk groups were presented in Figure S3A (all p < 0.001). Interestingly, the high-risk group had lower OS as compared to the low-risk group (Figure 6C, p = 0.002). The distribution of genomic subtypes, risk scoring groups, and OS status is shown in Figure 6D. The time-independent ROC curve further indicated the predicting accuracy of the prognostic model in the 1st, 2nd, and 3rd years (Figure 6E, all areas under the curve (AUCs) > 0.700). The AUC value was significantly higher as compared to 4 published signatures (Figure S2) (35–38). The distribution of risk score, OS status, and OS time is shown in Figure 6F. Notably, with the increasing risk score, the OS time of most HCC patients tended to be downregulated, which further verified the prognostic value of the model.

It was found that the relative scale of fraction for plasma cells, T-cell CD8, T-cell CD4 memory resting, T-cell CD4 memory activated, Tregs, macrophages M1, and DCs resting was obviously lower in the high-risk group than that in the low-risk group (Figure 7A, p < 0.01). By contrast, macrophages M0 and M2 were remarkably higher in the high-risk group than in the low-risk group (p < 0.01). In addition, all the ESTIMATE scores were significantly downregulated in the high-risk group (p < 0.001). We also analyzed the correlation between risk score and the relative levels of immune infiltration (Figure 7B). There were positive correlations between immune cells (B cells, CD4 cells, CD8 cells, DCs, macrophages, and neutrophils) and risk score. Further studies were performed to explore the relations between immune checkpoint expressions and risk scores (Figure 7C).

Based on the optimal TME value calculated by the mattools methods, all HCC patients were assigned to the high-TMB (n = 49) and the low-TMB groups (n = 301). The results of the K-M survival curve indicated that the OS rate for HCC patients with high TMB was significantly worse than that with low TMB (Figure 8A, p < 0.001). The prognostic value of risk score combined with TMB is shown in Figure 8B. We found that patients with low-TMB and low-risk scores had remarkably better prognostic status than others (p < 0.001). In addition, the oncoplots demonstrated the differences in mutated genes, mutation frequency, mutation types, and OS status between risk scoring groups (Figures 8C, D). Figure 8E indicates the results of the drug sensitivity test between risk score and six common chemotherapeutic drugs (including bosutinib, cyclopamine, dasatinib, gefitinib, metformin, and elesclomol) targeting HCC. Taken together, the risk score showed significant sensitivities to all chemotherapeutic drugs (p < 0.001). In the IMvigor210 cohort, the risk score of patients with stable disease or progressive disease (SD/PD) was remarkably lower than that of patients with complete response or partial response (CR/PR) (Figure S3B, p = 0.004).

The clinical prognostic value of this prognostic model was further validated in the FAHWMU cohort (n = 115). According to the above calculation formula, all HCC patients in the FAHWMU cohort were separated into the high-risk (n = 58) and low-risk groups (n = 57). The K-M survival curve indicated that the risk score had a good performance for predicting the OS outcomes (Figure 9A, p < 0.001). The time-independent ROC curve further demonstrated the prognostic value of the model (Figure 9B). It was found that the AUC value reached 0.816 in the 1st year, 0.922 in the 2nd year, and 0.928 in the 3rd year. The univariate Cox regression analysis of the risk score and clinical characteristics (including gender, age, TNM stage, tumor size, hepatitis B, Lymph node invasion, vascular invasion, perineural invasion, albumin, AFP, CEA, and CA199) was performed. Clearly, it was found that risk score, tumor size, AFP, and CEA were remarkably correlated with OS of HCC patients (Figure 9C, p < 0.05). Subsequently, they were used to construct a novel nomogram to better predict the prognosis of HCC (Figure 9D). In addition, the calibration curves were applied to validate the accuracy of the nomogram in the 1st, 2nd, and 3rd years (Figure 9E). The PH hypothesis analysis of the nomogram model is shown in Table S2.