Development and Validation of a Novel Prognosis Model Based on a Panel of Three Immunogenic Cell Death-Related Genes for Non-Cirrhotic Hepatocellular Carcinoma

Introduction

Globally, hepatocellular carcinoma (HCC) is the fifth most common cancer and the second leading cause of cancer-related deaths in men after lung cancer.¹ Most HCCs develop in patients with a history of cirrhosis. However, non-cirrhotic HCC (NCHCC) accounts for approximately 20% of all HCC cases.² Several studies have examined the diagnostic, therapeutic, and prognostic strategies for HCC.^3–6 Previous studies have demonstrated that the risk factors, pathogenetic backgrounds, clinical characteristics, tumor development, prognosis, and surveillance indicators vary between patients with cirrhotic HCC and those with NCHCC.^7–10 However, limited studies have examined the etiology, treatment, progression, and prognosis of NCHCC.^11,12 The diagnosis of NCHCC at an early stage is challenging owing to the lack of surveillance in non-cirrhotic patients, which may lead to delayed treatment.¹³ Currently, TNM staging is commonly regarded as a basis for formulating treatment strategies and predicting prognosis in NCHCC based on the macroscopic anatomical information of tumor tissues.¹⁴ However, TNM staging is inefficient with limited guiding significance as it does not account for molecular heterogeneity in tumors.¹⁵ Patients with the same TNM stage were reported to exhibit differential clinical outcomes.¹⁶ The landscape changes at the genomic level in tumor tissues can significantly affect the prognosis and therapy response of patients with cancer.^17–19 Therefore, there is a need to develop an early warning and prediction tool with molecular signatures to predict the prognosis of patients with cancer and supplement the TNM staging system. In the era of precision treatment, the prognostic factors of patients with NCHCC must be identified to facilitate clinical counseling and individualized prognosis prediction for NCHCC.

Immunogenic cell death (ICD) is a new type of cell death induced by various anticancer treatments, such as radiotherapy and chemotherapy. Additionally, ICD is a type of regulated cell death that occurs in infected and malignant cells and can activate both innate and adaptive immune responses by promoting the release of damage-associated molecular patterns (DAMPs).²⁰ The release of DAMPs can initiate a molecular cascade that induces the maturation of dendritic cells (DCs), which present antigens to cytotoxic T cells, activating tumor-specific CD8+ T lymphocytes to remove adjacent cells and promote immunological memory.^21,22 Additionally, DAMPs, such as CALR, TLR4, and HMGB1, which are immunogenic characteristics and biomarkers of ICD, are reported to be correlated with cancer prognosis.^23–26 ICD can improve the immunogenicity of tumors, promoting the transformation of cold tumors into hot tumors and potentiating the immunotherapeutic sensitivity of tumors.^27–29 Various cancer studies have used ICD-related markers to construct prognostic models for tumors, such as breast cancer and melanoma.^30,31 However, the functional role of ICD-related genes in NCHCC and their correlation with prognosis have not been elucidated. This study aimed to examine the functional role of ICDs in NCHCC and their ability to guide the therapeutic response of tumors and predict the prognosis.

This study screened NCHCC disease-related signature gene sets using weighted gene co-expression network analysis (WGCNA) with The Cancer Genome Atlas (TCGA) datasets, identified ICD-related molecular subtypes using the unsupervised clustering algorithm, and established a scoring model to predict the response to chemotherapy and immunotherapy. The ability of the model to predict the prognosis of patients with NCHCC was evaluated. A nomogram was constructed to rapidly evaluate the survival time. Additionally, the effects of GRAMD1C or TMC7 knockdown on the migration, invasion, and proliferation of HepG2 and Huh7 cell lines were examined using the wound healing, transwell, and Cell Counting Kit-8 (CCK8) assays, respectively. Furthermore, the correlation of the expression levels of GRAMD1C and TMC7 with survival in clinical patients was examined.

This study identified ICD-related genes associated with the prognosis of patients with NCHCC and constructed a nomogram to assess prognosis of NCHCC. This will aid clinicians to develop individualized counseling programs and treatments based on risk stratification and enhance the effectiveness of individualized treatment for patients with NCHCC.

Materials and Methods

NCHCC and Corresponding Adjacent Tissue Gene Expression Data Retrieval and ICD-Related Gene Identification

The clinicopathological information and raw RNA sequencing data were downloaded from TCGA database (https://portal.gdc.com (accessed on 10 June 2022)). Only the non-cirrhotic HCC group was selected based on the Ishak score (≤ 5).³² From previous studies, 57 ICD-related genes were rooted.^24,33

Identification of Differentially Expressed Genes (DEGs) and Enrichment Analysis

Data were normalized, and the DEGs between cancer and adjacent non-cancerous tissues were analyzed using the limma package in R software. DEGs were identified based on the following criteria: log₂ |fold change| > 0.15; adjusted p-value < 0.05. The Volcano plots were constructed to visualize the expression profiles of all DEGs. The top 50 DEGs were visualized using a heatmap. To analyze the underlying functions of potential targets, DEGs were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using the “ClusterProfiler” and “org.Hs.eg.db” R packages.

Identification of NCHCC Phenotype-Related Module Genes Using WGCNA

WGCNA is a data reduction and unsupervised classification method.³⁴ Based on the expression profile of NCHCC, the co-expression network was constructed using the “WGCNA” R package with the following parameters: corType = “Pearson”, mergeCutHeight = 0.25, minModuleSize = 60. Next, the correlation between the modules and phenotypes was illustrated using a heatmap. Modules with p < 0.05 were considered to be significantly correlated with NCHCC phenotype and included in subsequent analysis. The Venn diagram was used to show the intersecting genes among DEGs, NCHCC phenotype-related genes, and ICD-related genes.

Clustering Analysis of Samples

Based on the expression of these 5 ICD-related biomarkers, 135 NCHCC samples were typed using the ConsensusClusterPlus package. The effect of typing was assessed using principal component analysis. The optimal number of categories was determined according to the change in the area under the cumulative distribution function curve. The k value of the cluster category ranged from 2 to 9. DEGs between clusters were analyzed using the limma package.

Development of a Prognostic Model

All patients with NCHCC were equally randomized into a training set and a validation set. The training set data were subjected to univariable Cox regression analysis to identify the factors affecting the prognosis of patients. The significant factors identified in univariable Cox regression analysis were included in the stepwise variable multivariable Cox regression analysis to construct a risk score model for predicting disease prognosis. The prognostic risk score was calculated based on the regression coefficients (β) in the multivariable Cox regression model and the gene expression levels. A risk score was computed as follows: risk score = Σ_iCoefficient gene_i × expression of gene_i. The optimal cut-off points of risk score were obtained using X-tile software. Kaplan-Meier survival analysis revealed significant differences in the survival of different subgroups. To validate the accuracy of the model, the “time receiving operating characteristic (ROC)” package was used to calculate the area under the ROC curve (AUC) of patients with 1-year, 2-year, and 3-year overall survival (OS). The AUC values were directly proportional to the predictive performance.

Analysis of Intertumoral Immune Microenvironment

Single-sample gene set enrichment analysis (ssGSEA) algorithm was used to assess the immune function score of each patient using the “GSVA” R package. The ssGSEA enrichment score represents the relative abundance of each immune cell type or function. The unit distribution is normalized from 0 to 1. The expression of immune checkpoints in the two groups was examined. The correlation between risk score and immune checkpoints was evaluated using the “reshape2” R package. The CIBERSORT algorithm, which uses gene expression data to predict the proportions of 22 types of tumor-infiltrating immune cells (TIICs), such as T cells, B cells, macrophages, and natural killer (NK) cells, was used to estimate the differential immune cell infiltration status between the high-risk and low-risk NCHCC-TCGA groups using gene expression data. The RNAseq and maf files for all patients with NCHCC in TCGA cohort were obtained from cBioPortal. The “maftools” package was used to analyze and visualize genomic data. Additionally, the expression levels of HLA-related genes in different risk groups were analyzed.

Prediction of Response to Immunotherapy, Chemotherapy, and Targeted Therapy

The Genomics of Drug Sensitivity in Cancer database, which contains data on the sensitivity of 809 cell lines to 198 chemical compounds, was used to estimate the sensitivity of each patient to chemotherapy drugs. The half-maximal inhibitory concentration (IC50) is inversely proportional to the susceptibility to compounds. The Tumor Immune Dysfunction and Exclusion algorithm was used to evaluate the responses of patients in the two groups to immune checkpoint blockade therapies.

Clinical Value of the Prognostic Signature

To examine the clinical value of the risk score, univariable and multivariable Cox regression analyses were performed to estimate the independent effect of prognostic signatures (risk score, gender, grade, and stage) on the prognosis of TCGA-NCHCC cohort.

Construction of Nomograms

Nomograms were constructed to predict the prognosis of patients with NCHCC based on age, stage, and risk score. The ROC and calibration curves were plotted using the “timeROC” and “rmda” R package, respectively, to assess the diagnostic value of the nomogram.

Immunohistochemical (IHC) and Tissue Microarray (TMA) Analyses

The absence of liver cirrhosis was confirmed after a complete pathological examination of the resected specimen. TMAs comprised 53 NCHCC samples from different cases. Sample collection was approved by the Research Ethics Committee of The First Affiliated Hospital of Chongqing Medical University (Reference number: 2022–K518). The sections were warmed, deparaffinized, rehydrated, subjected to heat-induced antigen retrieval, blocked, and washed. Immunoreactive signals were detected using an Envision/horseradish peroxidase system (Dako-Cytomation, Glostrup, Denmark). All sections were counterstained with hematoxylin and evaluated under a light microscope.

Cell Culture and Transfection

HepG2 and Huh7 cells were purchased from MeisenCTCC (Meisen Chinese Tissue Culture Collections) and cultured in high-glucose Dulbecco’s modified Eagle’s medium (DMEM) (Meilunbio, Dalian, China) supplemented with 10% fetal bovine serum (FBS) (Meiluncell, Dalian, China) and 1% penicillin-streptomycin mixture at 37°C and 5% CO₂. Cells were transfected with 50 nM of specific short-interfering RNA (siRNA) or a negative control siRNA using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) in opti-MEM (Meiluncell, Dalian, China). The cells were cultured in high-glucose DMEM supplemented with 10% FBS after transfection.

RNA Extraction and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Total RNA was extracted from cells using Trizol reagent (Takara, Beijing, China) and reverse-transcribed into complementary DNA using the PrimeScript^TM RT reagent kit with genomic DNA eraser (Takara, Beijing, China), following the manufacturer’s instructions. qRT-PCR analysis was performed with SYBR Green qPCR Master Mix (MedChemExpress, Monmouth Junction, NJ, USA). GAPDH served as the internal reference gene. Various primers were used to perform qRT-PCR analysis. The mRNA levels were determined using the 2^−∆∆CT method. Samples were analyzed in triplicates.

Western Blot Analysis

Total proteins were extracted from cells using the radioimmunoprecipitation assay protein extraction reagent (Meilunbio, Dalian, China) at 48 h post-transfection. The proteins were quantified using the bicinchoninic acid protein assay kit (Beyotime, Jiangsu, China). Proteins were subjected to sodium dodecyl sulfate-polyacrylamide gel electrophoresis using a 7.5% gel. The resolved proteins were transferred to a polyvinylidene difluoride membrane (0.22µm, Merck Millipore, Germany). The membrane was blocked, incubated with primary antibodies, washed, and visualized using an enhanced chemiluminescence detection system (Bio-Rad Laboratories). ACTB was used as an internal control. Western blot analyses were performed in triplicate.

Transwell and Wound Healing Assays

Wound healing and transwell assays were performed on day 2 post-transfection. Cells were seeded in a 6-well plate and cultured until they reached 100% confluency and formed a confluent monolayer. A sterile p200 pipette tip was used to scratch the monolayer. Cellular debris was removed by gently washing the monolayer with phosphate-buffered saline. Points near the scratch were marked to serve as reference points. Wound closure was examined under a microscope at 0 and 24 h. The free wound area was measured using the ImageJ software. The migration ability of cells was determined based on the closure of the wound gap.

The polycarbonate filter membrane at the bottom of the transwell chamber was covered with 70 µL Matrigel (Solarbio, Beijing, China), air-dried, and maintained overnight in a laminar flow cabinet. At 24–48 h post-transfection, HepG2 and Huh7 cells were seeded in the upper transwell chamber (pore size, 12 μM; diameter, 8 mm; BIOFIL, Guangzhou, China) at a density of 4×10⁴ cells per chamber. DMEM supplemented with 30% FBS was added to the lower compartment (600 µL per well). After incubation at 37°C and 5% CO₂ for 48 h, cells that failed to migrate were gently wiped off. The membrane was fixed with paraformaldehyde for 15 min and stained with 0.5% crystal violet. The number of invading cells in five random fields was counted under an inverted microscope (200×) (OLYMPUS CKX53, Japan). All experiments were performed in triplicate.

CCK8 Assay

HepG2 and Huh7 cells in the logarithmic growth phase were seeded in a 96-well plate (5 ×10³ cells/well) at 8 h post-transfection and cultured for 24, 48, and 72 h. Next, the cells were incubated with 10 μL CCK8 solution (Beyotime Biotechnology, Shanghai, China) for 1 h and processed using the CCK8 kit, following the manufacturer’s instructions. The absorbance of the reaction mixture at 450 nm was recorded using Varioskan Flash (Thermo Scientific, USA).

Statistical Analyses

Bioinformatics analyses were performed using the R project (3.6.3). Risk comparison between unpaired samples was performed using the Mann–Whitney U-test. Spearman correlation test was used for various analyses, such as ssGSVA and CIBERSORT analysis. The hazard ratio and 95% confidence interval were computed using Cox regression analyses. Kaplan-Meier survival curves were compared using the Log rank test. All statistical analyses were performed using the SPSS22.0 statistical software (SPSS Inc., Chicago, IL, USA). Data are expressed as mean ± variance. Means between two groups were compared using the t-test, whereas those between multiple groups were compared using one-way analysis of variance. Enumeration data were compared using the χ² test. Differences were considered significant at P < 0.05.

Results

Analysis of DEGs in NCHCC

The study workflow is shown in Figure 1. The data of 377 patients with HCC, including RNA sequencing expression profiles and corresponding clinical information, were downloaded from TCGA. Of these, 159 cases were excluded due to the lack of Ishak score information (Figure S1). The remaining 218 cases were screened based on their Ishak scores. The patients with scores ≤ 5 were included in the non-cirrhotic group (n = 133).³⁵ Finally, 133 non-cirrhotic HCC and 29 corresponding adjacent tissue data were examined. The characteristics of the patients and lesions are summarized in Table 1. In this study, 3167 DEGs between cancer and adjacent tissues were identified (1809 upregulated and 1358 downregulated genes) (Figure 2a). Functional enrichment analyses were performed to determine the biological classification, function, and pathway of DEGs. GO analysis (Figure S2a) revealed that in the biological processes category, DEGs were enriched in response to the proteasomal protein catabolic process, mitotic cell cycle phase transition, negative regulation of phosphate metabolic process, negative regulation of phosphorus metabolic process, and regulation of protein catabolic process. Meanwhile, in the cellular component category, DEGs were enriched in cell-substrate junction, nuclear speck, chromosomal region, chromosome, centromeric region, and outer membrane. In the molecular function category, DEGs were enriched in kinase regulator activity, protein kinase regulator oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, iron ion binding, and ribonuclease activity. KEGG pathway analysis (Figure S2b) revealed that DEGs were enriched in cellular senescence, viral carcinogenesis, alcoholism, cell cycle, and TNF signaling pathway.

Table 1 Characteristics of Patients with Non-Cirrhotic Hepatocellular Carcinoma in the Cancer Genome Atlas Database

Figure 1 The flowchart of the study.

Figure 2 Identification of genes associated with ICD and the phenotypic and expression profiles of NHCC. (a) The heatmap of the top 50 differentially expressed genes between the two groups. (b–d) Identification of key gene modules using weighted gene co-expression network analysis. (b) Sample clustering chart. (c) Module clustering dendrogram. (d) Heatmap of the relationship between gene modules and traits (the non-cirrhosis type was used as the phenotype). (e) The Venn diagram of DEGs, ICD-related genes, and key module genes. Abbreviations: ICD, immunogenic cell death; NCHCC, non-cirrhotic hepatocellular carcinoma; DEGs, differently expression genes.

Screening of modules corresponding to NCHCC phenotype and identification of intersection genes to identify non-cirrhosis-related genes in NCHCC, TCGA-NCHCC dataset was subjected to WGCNA (Figure 2b–d). Genes corresponding to brown, pink, and gray modules were selected for subsequent analysis. After the identification of cancer-related, NCHCC phenotype-related, and ICD-related genes, 5 intersection genes were obtained from the Venn diagram (Figure 2e).

Acquisition of Intersection Genes

ICD-related gene-based clustering analysis, molecular typing, and prognostic assessments consensus clustering analysis was performed with 135 NCHCC cases based on the expression profiles of 5 intersection genes, which revealed the following two subtypes: subtype A (65 cases) and subtype B (68 cases) (Figure 3a). Kaplan-Meier curves revealed that the OS of patients in subtype A was higher than that of patients in subtype B, suggesting that subtype A was associated with improved prognosis (Log rank test, p < 0.001; Figure 3b). Next, 186 DEGs between subtypes A and B were identified using the limma package and the cut-off criteria (Figure 3c). In total, 147 prognosis-related DEGs were screened from the intersection of subtype-related and tissue-related DEGs between cancer and adjacent tissues using the Venn diagram (Figure 3d).

Figure 3 ICD-related-gene-based clustering analysis, molecular typing and prognostic assessments. (a) Consensus matrix based on the expression of five genes in the NCHCC cohort at k = 2. (b) K-M survival analysis between two subtypes. (c) The heatmap of the top 50 DEGs between the two subtypes. (d) The Venn diagram of subtypes DEGs and NCHCC DEGs. Abbreviations: ICD, immunogenic cell death; NCHCC, non-cirrhotic hepatocellular carcinoma; DEGs, differentially expressed genes; K-M, Kaplan-Meir.

Construction and Validation of the Prognosis Model

Univariable and multivariable Cox regression analyses were performed to identify the prognostic factors for OS based on 147 prognosis-related DEGs. Univariable Cox regression analysis revealed 18 prognosis-related factors in NCHCC (Figure 4a). Next, these 18 factors were subjected to stepwise variable multivariable Cox regression analysis to develop a risk score model for predicting prognosis. The panel comprising TMC7, GRAMD1C, and GNPDA1 was a significant predictor of OS in patients with NCHCC. The risk score was calculated as follows: risk score = (1.04 × expression level of TMC7) + (−1.49 × expression level of GRAMD1C) + (2.15 × expression level of GNPDA1). Based on the risk score and optimal cut-off points, the training set was divided into high-risk and low-risk groups. Kaplan-Meier curves revealed that the OS of patients in the low-risk group was higher than that of patients in the high-risk group (Figure 4b and c). The model exhibited good performance in distinguishing the low-risk group from the high-risk group. The AUC values for 1-year, 2-year, and 3-year survival were 0.805, 0.666, and 0.674, respectively. As shown in Figure 4d–g, the performance of the risk score model was good in both the validation set and total TCGA-NCHCC samples.

Figure 4 Construction and validation of the prognosis model. (a) Forest plot of univariable Cox regression results for the 18 DEGs. Red represents HR>1, green represents HR<1. (b and c) Kaplan-Meier survival curves of the high-risk and low-risk groups in the training set. The survival ROC curve of the training set. (d and e) Kaplan-Meier survival curves for the high-risk and low-risk groups in the validation set. The survival ROC curve of the validation set. (f and g) Kaplan-Meier survival curves of the high-risk and low-risk groups in all sets. The survival ROC curve of all sets. Abbreviations: DEGs, differentially expressed genes; HR, hazard ratio; ROC, receiving operating characteristic.

Analyses of Tumor Microenvironment in Different Risk Groups

To understand the heterogenicity of the tumor microenvironment of the high-risk and low-risk groups, immune-related analyses, including ssGSEA and CIBERSORT, were performed. ssGSEA revealed that several signaling pathways, such as APC_co_stimulation, check−point, MHC_class_I, and type_II_IFN_reponse were mainly enriched in the high-risk group (Figures 5a and S3a). The CIBERSORT deconvolution approach was used to evaluate the relative abundance of 22 TIICs in the high-risk and low-risk groups. B cells memory, T cells regulatory (Tregs), Dendritic cells resting, Mast cells resting were markedly enriched in the high-score group (Figure 5b). Next, the TIICs enriched with each signature gene in the risk score model were analyzed. The expression of GNPDA1 was positively correlated with the proportion of Mast cells resting. Meanwhile, the expression of TMC7 was positively correlated with the proportion of Tregs and Dendritic cells resting but negatively correlated with the proportion of NK cells activated (Figure 5f). The proportion of Mast cell resting was moderately correlated with the risk score (Figure S3b). Next, the correlation between immune checkpoint expression and risk score was examined. The expression levels of immune checkpoints were significantly upregulated in the high-risk group (Figures 5c and S3c). Additionally, the expression levels of 24 immune checkpoint genes (such as TNFRSF18, NRP1, and CD44) were positively correlated with the risk score.

Figure 5 Analyses of tumor microenvironment in different risk groups. (a) The fraction of signaling pathways in the high-risk and low-risk groups. (b) The fraction of tumor-infiltrating immune cells in the high-risk and low-risk groups. (c) The expression levels of immune checkpoints in the high-risk and low-risk groups. (d) The waterfall graphs of the most frequently mutated genes in different risk groups. (e) The expression levels of HLA-related genes in the high-risk and low-risk groups. (f) The association between tumor-infiltrating immune cells and biomarker genes. *P < 0.05, **P < 0.01, and ***P < 0.001.

Tumor mutational burden (TMB) is closely associated with cancer immunotherapy.³⁶ The most common alteration in the high-risk score group was CTNNB1 mutations (42.9%), followed by TP53 (32.1%) and MUC16 and TTN (>20%) mutations. In the low-risk score group, the most common alteration was TP53 mutation (31.7%), followed by TTN (25%) and CTNNB1 (20%) mutations (Figure 5d). The expression levels of HLA-related genes, such as HLA-A, HLA-DMB, HLA-DMA, HLA-DOA, HLA-DPA1, HLA-DPB2, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DQB1, HLA-DQB2, HLA-DRA, HLA-DRB1, and HLA-DRB6 varied between different groups (Figures 5e and S3d).

Prediction of the Drug Sensitivity of NCHCC

The IC50 values of common chemotherapy drugs, small-molecule drugs, or targeted drugs against HCC were analyzed in the high-risk and low-risk groups.^37–40 The IC50 values of vinblastine, axitinib, GDC-0449, methotrexate, and docetaxel against tumors in the high-risk group were significantly higher than those against tumors in the low-risk group (Figure 6a–e, p < 0.05). In contrast, the IC50 values of bleomycin, cisplatin, doxorubicin, epothilone B, gemcitabine, mitomycin C, obatoclax, mesylate, and sorafenib against tumors in the high-risk group were lower than those against tumors in the low-risk group, indicating that these prognostic-related genes may influence the therapeutic efficacy of drugs against HCC (Figure 6f–m, p < 0.05).

Figure 6 (a–m) IC50 values of common chemotherapy drugs, small molecule drugs, or targeted drugs for HCC in different risk groups. Abbreviations: IC50, half-maximal inhibitory concentration; HCC, hepatocellular carcinoma.

Clinical Value of the Prognostic Signature

The evaluation of the independent prediction ability of the risk signature using univariable and multivariable Cox regression analyses of signature and other common prognostic factors revealed that the risk score and tumor grade were associated with the prognosis of NCHCC in TCGA-NCHCC cohort (Figure 7a and b). Analysis of the distribution of the risk score in different subgroups revealed that the risk score was directly proportional to the grade stage (Figure S3e–f). The Kaplan-Meier curves of the subgroups revealed that the prognosis of the high-risk score group was poor when compared with that of the low-risk score group in the early and middle stages of the disease (Figure 7c–g). Similar patterns were observed with M0, T1, and N0 tumors. These findings indicate that patients with NCHCC of the same clinical stage may exhibit differential prognosis. The risk score model developed in this study should be used to distinguish patients with poor prognoses to optimize treatment plans and consequently achieve an improved prognosis.

Figure 7 Clinical value of the prognostic signature. (a) Forest plot of univariable Cox regression results for 5 prognostic factors. (b) Forest plot of multivariable Cox regression results for 5 prognostic factors. (c–g) Kaplan-Meier survival curves of the high-risk and low-risk groups according to the subgroups.

To predict the probability of 1-year and 3-year survival, a nomogram was constructed with age, stage, and risk score (Figure 8a). ROC curves were used to evaluate the diagnostic value of the nomogram (Figure 8b). The AUC values of the nomogram for 1-year and 3-year survival were 0.799 and 0.737, respectively. These results along with the calibration curves demonstrated the potent diagnostic value of the nomogram (Figure 8c).

Figure 8 Prognosis prediction model for patients with NCHCC. (a) Nomogram for 1-year and 3-year OS rates. (b) Time-dependent ROC curves of the nomogram. (c) Calibration plots of the nomogram. *P < 0.05; ***P < 0.001. Abbreviations: NHCC, non-cirrhotic hepatocellular carcinoma; OS, overall survival; ROC, receiver operating characteristic.

Effects of TMC7 and GRAMD1C on the Malignant Biological Behavior of Liver Cancer Cells

A previous study demonstrated that GNPDA1 can promote the proliferation, migration, and invasion of HCC cells and inhibit apoptosis.⁴¹ To evaluate the role of TMC7 and GRAMD1C in the migration, invasion, and proliferation of liver cancer cells, TMC7 and GRAMD1C were knocked down in HepG2 and Huh7 cell lines using siRNA. The knockdown efficiencies of siRNAs were confirmed using qRT-PCR and Western blotting analyses (Figure 9a–c). The effect of TMC7 and GRAMD1C knockdown on HCC cell malignant behavior was examined using the CCK8, transwell, and wound healing assays. As shown in Figure 9d and e, the migration of TMC7 knockdown HCC cells was significantly downregulated, whereas that of GRAMD1C knockdown HCC cells was significantly upregulated. The results of the transwell assay revealed that tumor cell invasion was significantly suppressed upon TMC7 knockdown but significantly enhanced upon GRAMD1C knockdown (Figure 9f and g). However, the results of the CCK8 assay revealed that the knockdown of TMC7 and GRAMD1C did not affect the proliferation of HepG2 and Huh7 cells (Figure 9h and i).

Figure 9 TMC7 and GRAMD1C play a key role in cell migration and invasion abilities of HepG2 and Huh7 cells. (a) qRT-PCR analysis revealed that the TMC7 and GRAMD1C mRNA levels in the si-TMC7-transfected and si-GRAMD1C-transfected groups, respectively, were significantly downregulated when compared with those in the si-NC-transfected group. (b and c) Western blotting analysis demonstrated that the TMC7 and GRAMD1C protein levels in the si-TMC7-transfected and si-GRAMD1C-transfected groups, respectively, were significantly downregulated when compared with those in the si-NC group. (d and e) Effect of TMC7 and GRAMD1C knockdown on the migration ability of HepG2 and Huh7 was evaluated using the wound healing assays. (f and g) Effect of TMC7 and GRAMD1C knockdown on the invasion ability of HepG2 and Huh7. (h and i) Effect of TMC7 and GRAMD1C knockdown on the proliferation ability of HepG2 (h) and Huh7 (i) cells was evaluated using the CCK-8 assay. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. Abbreviations: NC, negative control; si, small interfering RNA; qRT-PCR, quantitative real-time polymerase chain reaction; CCK-8, cell counting kit-8.

The expression levels of TMC7 and GRAMD1C in NCHCC tissues and corresponding tumor-adjacent tissues from the TMA were examined using IHC analysis. Tumor tissues exhibited upregulated TMC7 expression and downregulated GRAMD1C expression (Figure 10a and b). The upregulated expression of TMC7 and the downregulated expression of GRAMD1C predicted poor outcomes in patients with NCHCC (Figure 10c and d).

Figure 10 Immunohistochemical analysis of two genes in NCHCC tissues and corresponding tumor-adjacent tissues and the Kaplan-Meier curves of patients according to the expression of these two genes. (a and b) Immunohistochemical analysis of two genes in NCHCC tissues and corresponding tumor-adjacent tissues in TMA. (c and d) Kaplan-Meier analysis of overall survival of patients according to the expression of these two genes. Abbreviations: NCHCC, non-cirrhotic hepatocellular carcinoma; TMA, tissue microarray.

Discussion

In the early stages of ICD, the “do not eat me” and “eat me” signals are mutually restraining. Additionally, tumor cells can be effectively recognized and phagocytosed by DCs and macrophages, triggering antitumor immune responses.⁴² Therefore, ICD can be an effective pathway to activate the immune response against cancer and can determine the long‐term success of anticancer therapies. In the direct-acting antivirals era, the incidence of non-cirrhotic HCC is increasing, whereas that of cirrhotic HCC is decreasing.⁴³ NCHCC is not diagnosed at early stages as it does not exhibit symptoms during early stages. The pathogenesis of NCHCC is distinct from that of cirrhotic HCC. Thus, HCC therapy may not be the optimal treatment for NCHCC.¹² The activation of ICD can improve the immunogenicity of tumors and enhance the immunotherapeutic sensitivity of tumors. The correlation between ICD-related genes and NCHCC prognosis has not been previously examined. Therefore, this study independently examined the prognosis of NCHCC and the prediction of drug therapy based on ICD-related genes.

In this study, the correlation between ICD and NCHCC prognosis was evaluated based on RNA sequencing data from TCGA. NCHCC cases were classified into two subgroups based on key genes associated with ICD, non-cirrhotic phenotype, and cancer tissue using consensus clustering. The prognosis of the two subgroups significantly varied, indicating the presence of two distinct molecular subtypes in NCHCC. A risk score model was constructed for the prediction of NCHCC prognosis using a panel comprising TMC7, GRAMD1C, and GNPDA1. Compared with several other immune-related prognostic models comprising multiple parameters, the risk score model developed in this study comprised 3 independent prognostic genes, which can be clinically useful.

GNPDA1 can catalyze the conversion of glucosamine 6-phosphate to fructose 6-phosphate and consequently provide raw materials for glycolysis, promoting cancer progression.⁴⁴ Additionally, GNPDA1 plays important roles in cell proliferation, migration, and invasion.^41,45,46 The upregulation of GNPDA1 is associated with poor prognosis in patients with HCC.⁴¹ TMC7 is a member of a gene family predicted to encode transmembrane proteins.⁴⁷ The expression of TMC7 is upregulated in pancreatic carcinoma, rectal cancer, and head and neck squamous cell carcinoma, contributing to tumor progression and metastasis.^48–50 GRAMD1C, a negative regulator of starvation-induced autophagy, can regulate autophagy initiation and mitochondrial bioenergetics.⁵¹ The downregulation of GRAMD1C was associated with advanced clinicopathological characteristics (advanced clinical stage, distant metastasis, and advanced histological grade), poor prognosis, and decreased OS in kidney renal clear cell carcinoma.⁵² In this study, GNPDA1 and TMC7 were determined to be adverse prognostic factors for NCHCC, whereas GRAMD1C was a favorable prognostic factor, which is consistent with the results of previous studies.

HCC cells reshape the tumor microenvironment through various mechanisms, enabling tumor cells to evade immune surveillance by decreasing the number of T cells and the activity of NK cells, resulting in the upregulation of tumor cell proliferation and metastasis.⁵³ The poor prognosis of patients with HCC was correlated with the increased proportions of Tregs.⁵⁴ Analysis of immune cell infiltration in cancer tissues revealed that the expression of TMC7 was positively correlated with the proportion of Tregs but negatively correlated with the proportion of activated NK cells. This suggested that the poor prognosis of the high-risk group may be related to immunosuppression.

Patients with NCHCC were classified into high-risk and low-risk groups using the risk score model. The immune microenvironment was significantly different between the groups. A previous study demonstrated that type II IFN (γ-IFN) can mediate immune escape and immunosuppression to promote tumor progression.⁵⁵ Compared with those in the low-risk group, the immune checkpoint gene expression levels and the enrichment of immune pathways were significantly upregulated and the levels of type II IFN response were significantly downregulated in the high-risk group. This can explain the increased sensitivity of the high-risk group to small-molecule and immunotherapy drugs.

The HLA system mainly comprises the following two classes: major histocompatibility complex (MHC) class I, comprising HLA-A, HLA-B, and HLA-C; MHC class II, comprising HLA-DP, HLA-DQ, and HLA-DR. MHC class I, which is expressed on the surface of all nucleated cell types, is recognized by CD8+ T cells.^56,57 Meanwhile, MHC class II is expressed mainly in antigen-presenting cells, including DCs and macrophages, which present exogenous antigens. A previous study reported that the upregulation of MHC class II suppresses tumor-infiltrating lymphocytes and tumor immune response.⁵⁸ In this study, analysis of HLA-related genes indicated that the expression levels of several MHC class II genes (HLA-DPA1, HLA-DPB2, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DQB1, HLA-DQB2, HLA-DRA, HLA-DRB1, and HLA-DRB6) in the high-risk group were significantly higher than those in the low-risk group. We speculated that the high-risk group may be associated with an immunosuppressive state, which is correlated with poor prognosis. TMB is a predictive biomarker for immunotherapy in multiple cancers.⁵⁹ In this study, TMB analysis revealed that the mutation frequency of CTNNB1 in the high-risk group was significantly higher than that in the low-risk group. However, the differences in the frequency of other mutations were minimal between the two groups. HCC with a high proportion of CTNNB1 mutation, which is related to the inhibition of immune cells and progression of tumors, is highly likely to be insensitive to immunotherapy.^60,61 Drug sensitivity analysis revealed the sensitivity of tumors in the high-risk group to multiple anticancer drugs was distinct from that of tumors in the low-risk groups. Consequently, appropriate anticancer drugs can be selected based on the risk score for patients with NCHCC.

Clinicians often classify and grade patients with tumors according to tumor size, the degree of invasion, lymph node involvement, and distant metastasis (TNM) to decide the treatment strategy.⁶² However, this study demonstrated that even patients with the same stage tumor may exhibit differential prognostic outcomes. Therefore, patients with the same stage should be further stratified according to the risk score model developed in this study and subjected to aggressive treatment. For early-stage disease, curative intent treatments include surgery, radiofrequency ablation, and liver transplant.^63,64 Meanwhile, for advanced-stage disease, conventional cytotoxic chemotherapy is ineffective with limited clinical benefits. Treatment with multikinase inhibitors is expected to enhance the OS of patients with advanced HCC by only a few months.^6,64 Accurate risk assessment in the early stage of liver cancer will aid in improving therapeutic response and prognosis. The risk score model developed in this study can assess the risk level in NCHCC at early stages (stage 1, M0, T1, and N0). This strategy will enable clinicians to accurately predict the prognosis of patients, devise treatment plans, and communicate with patients and their families about the treatment plan rather than only evaluating the clinical condition of patients based on rough indicators (such as the TNM stage).

The predictive nomogram integrates the molecular prediction model at the micro level and the staging index of the TNM stage at the macro level and provides a convenient and precise method for clinicians to evaluate the 1-year and 3-year prognosis of patients with NCHCC. The mortality rate and corresponding survival rate within a certain period can be evaluated by scoring three prognostic indicators.

This study further verified the correlation between the expression of variables in the risk model and the malignant biological behavior of HCC cells in vitro. The results of in vitro experiments were consistent with those of bioinformatics analysis. The mechanism of key genes affecting the migration and invasion abilities of liver cancer cells will be analyzed in future studies.

This study has some limitations. First, as the complete Ishak score information was lacking in other public databases, only internal verification was performed based on TCGA datasets. More cases must be used for external validation of the signature genes developed in this study. Second, in histopathological verification of signature gene expression, most of the included cases had a history of hepatitis B, which can be attributed to the high prevalence of hepatitis B in China.⁶⁵ However, other types of cases will be examined in the future to validate the conclusions of this study. Third, the model developed in this study can be used to predict patient survival. However, the mechanisms through which gene interactions in the model affect tumor progression have not been investigated. We plan to expand our research in metabolomics, proteomics and genomics to provide insight into downstream mechanisms.

In conclusion, bioinformatics analyses were performed to generate a three-gene predictive risk score model to predict the outcomes of patients with NCHCC. A novel nomogram was developed to predict the survival outcomes of these patients. Risk stratification will enable the administration of optimal therapy for patients with NCHCC.