A novel immune-associated prognostic signature based on the immune cell infiltration analysis for hepatocellular carcinoma

Data collection and grouping

Transcriptional group RNA-sequencing information and clinical messages were obtained from TCGA. The dataset comprised data on 371 HCC and 50 non-tumor tissues. The LIHC transcriptome data was grouped from TCGA based on ICI using the ssGSEA. This analysis was used to assess the infiltration levels of various immune cells, immune-associated signaling pathways, and functional activities. According to the median of ssGSEA results, the HCC samples were split into high-ICI (IMM_H) and low-ICI (IMM_L) groups. To determine the classification, the ESTIMATE algorithm was used to calculate various scores associated with the tumor microenvironment. These scores included Tumor Purity, Immune, ESTIMATE, and Stromal Scores. The box chart examination indicated that IMM_H was positively associated with these scores. No ethical approval was needed in this study since the data were retrieved from a public database.

Effectiveness of immune grouping

The ESTIMATE algorithm was used to identify differentially expressed genes (DEGs) in the HCC data. The algorithm determines the interstitial, estimated, immune, and tumor purity scores. The results were visualized by generating clustering heatmaps and statistical heatmaps (Supplementary Figures 1 and 2). We utilized human leukocyte antigen (HLA) genes to investigate the differences between the two types of immune cells. Additionally, the CIBERSort was used to measure the immune cell types and their association with the DEGs. These analyses comprehensively described the differences between the two groups in terms of ICI and gene expression levels.

Differentially expressed genes

After dividing the HCC samples into two groups, R software was used to identify DEGs related to immunity. A two-way Venn analysis was used to identify overlapping DEGs between the two groups. Functional enrichment of the DEGs was performed using Gene Ontology (GO) [20] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [21] using the clusterProfiler package [22] in R with p<0.05 indicating statistical significance, and thus suggesting a potential biological value of the enriched terms in the DEGs.

PPI and immune-related gene

The protein-protein interaction (PPI) network was constructed with the STRING database [23] and visualized with Cytoscape software (Version 3.7.0, National Resource for Network Biology, USA). Only those nodes with interaction confidence intervals that were more significant than 0.7 were included [24]. Thus, the interactions between IRGs and potential functional relationships could be determined. Next, the survival package in R was used to evaluate IRGs in the training dataset using univariate Cox regression analysis, to determine the association between IRG expression and patient survival. We further screened the IRGs using a log-rank test, identifying genes significantly impacting patient prognosis. Then, a multi-factor Cox proportional hazards regression (PHR) model was developed to determine the risk score for patient survival in terms of the IRGs.

Immune-related gene prognostic signature

The clinical HCC data from TCGA were included in the univariate Cox PHR evaluation. The significance was set at p<0.05 to identify IRGs significantly connected with patient survival. Next, multivariate Cox regression was used to identify prognostic markers (age, sex, tumor grade, tumor status, tumor stage, and risk score) using the identified IRGs. The risk measures were determined for each sample depending on the outcomes of the above analysis. A Kaplan–Meier (KM) survival analysis was performed to equate the distinctions between the high and low-risk categories based on survivability to predict the risk scores and prognosis. Then, the risk score was examined to determine gene expression levels using the formula:

External validation

The International Cancer Society database was used as the validation dataset to confirm the predictive value of the immune-related gene risk score signature. According to the median risk measurement, the HCC patients were categorized into high- and low-risk categories in the validation cohort. A KM survival analysis was conducted in the validation cohort. Additionally, receiver operating characteristic (ROC) curves facilitated the assessment of risk measurement presentation signature in predicting patient outcomes.

Correlation analysis of TICs

The infiltration data from various cell types in the immune system was obtained from the Tumor Immune Assessment Resource (TIMER) databank. Pearson correlation analysis was also used.

Statistical analysis

R software (Version 3.6.1, The R Foundation, initially written by “R & R” of the Statistics Department of the University of Auckland.) was used for all analyses, with p-value <0.05 or |log₂FC|>1 for assessment of significance.

Analysis process

The analysis had the following steps (Figure 1):

• Sample collection: Data on a total of 371 HCC and 50 adjacent tissues were retrieved from TCGA.

• ICI analysis: The HCC tissues were grouped based on ICI levels: IMM_H (n=154) and IMM_L (n=217).

• Identifying differentially expressed IRGs: A total of 177 IRGs were identified based on the screening criteria using R.

• Functional enrichment analysis: GO and KEGG evaluations determined the probable characteristic mechanisms and signals associated with identified IRGs.

• Protein-protein interaction (PPI) network construction: Overlapping DEGs were used to build a PPI network.

• Cox regression analysis and prognosis construction of risk signature: This was used to verify the genes linked with prognosis using the PPI network DEGs. This led to the construction of a prognosis risk signature.

• Assessment of the prognostic signature: This was evaluated using survival analysis, ROC analysis, principal component analysis (PCA), GSEA, and TIC correlation analysis.

• Validation of the signature: The prognostic signature was validated with an external database.

Figure 1:

Workflow of the study.

Construction and verification of HCC groupings

An unsupervised hierarchical clustering algorithm was used to divide the 371 HCC cases retrieved from TCGA into IMM_H (n=154) and IMM_L (n=217) according to the median. ssGSEA was used to evaluate 24 types of infiltrating immune cells (Figure 2A). Conversely, IMM_L showed a positive correlation with tumor purity (Figure 2B). Thus, the IMM_H group had higher immunity features and lower tumor purity than the IMM_L group (p<0.05). Furthermore, the HLA family expression was higher in the IMM_H group than in the IMM_L group (Figure 2C). The CIBERSort algorithm was used to verify ICI in the two groups. The IMM_H group showed a greater variety of immunogenic cells (Figure 2D). The IMM_H group had higher ICI, lower tumor purity, and enhanced IRG expression, suggesting a more favorable immune microenvironment for novel immunotherapy interventions.

Figure 2:

Development and verification of HCC grouping. (A) Grouping information of immune cells. (B) Between the groups, the box plot indicated a substantial distinction in Tumor Purity, ESTIMATE, Immune, and Stromal scores (p<0.01). (C) The levels of HLA homeostatic genes in IMM_H (red) were significantly greater than in IMM_L (blue) (p<0.01). (D) CIBERsort analysis showing the percentage distinction between the IMM_H (red) and the IMM_L (blue) for each immune cell.

Immune-related genes

We conducted a differential expression analysis between samples acquired from TCGA. We identified 1,833 DEGs, including 1,218 upregulated and 615 downregulated genes (Figure 3A) using the above criteria. Furthermore, a similar analysis was conducted between the two groups in the HCC samples. A total of 484 IRGs were identified with the same standards, including 398 upregulated and 86 downregulated genes (Figure 3B). A two-way Venn analysis found 177 common genes between the tumor and normal groups and the IMM_H and IMM_L groups (Figure 3C). These 177 p genes demonstrated the DEG intersection in the tumor and ICI analyses. Enrichment analyses were used to assess the molecular processes and pathways associated with the IRGs. GO analysis indicated that the IRGs were primarily enriched in the biological process category. In terms of cellular components, the IRGs were enriched in the extracellular matrix, collagen-containing extracellular matrices, and collagen trimer categories. In the molecular function category, enrichment was seen in the extracellular matrix structural constituents, receptor–ligand interactions, and regulatory activities (Figure 3D). Additionally, KEGG analysis showed that the IRGs were primarily associated with receptor interaction of cytokines, digestion, and absorption of proteins, and chemokine signaling (Figure 3E). A set of IRGs were found to be differentially expressed in HCC tissues. This information enhances the understanding of the immunological characteristics of HCC, generating novel insights into precision drug targets and prognostic markers.

Figure 3:

Analysis of IRGs. (A) Volcano plot indicating 1,218 upregulated and 615 downregulated DEGs between HCC and normal tissues. (B) Identification of 398 upregulated and 86 downregulated genes between IMM_H and IMM_L of HCC; red dots indicate upregulation and green dots downregulation (fold change>1, p=0.05). (C) 177 IRGs were identified by R software. (D and E) GO and KEGG enrichment of 177 IRGs, with p<0.05 used to show significant enrichment.

Immune-related gene prognostic signature

The STRING database and Cytoscape software were used to construct the PPI network to investigate the potential mechanism of the 177 DEGs. From the network, 82 key DEGs were identified with probable significant roles in HCC (Figure 4A). Univariate Cox regression analysis was conducted to identify 11 prognosis-linked IRGs (Figure 4B). Subsequently, multivariate Cox regression analysis of these 11 IRGs found that four IRGs, namely, C7, CTSV, MMP1, and VCAN, could be independent predictors of HCC patient outcomes (Figure 4C).

Figure 4:

PPI network and Cox regression analyses. (A) The interaction network was constructed using nodes with an interaction confidence of >0.70. (B) The 82 IRGs underwent univariate Cox regression analysis, using p<0.05 as the level of significance. (C) Multivariate Cox regression analysis for determination of prognosis-associated genes.

We measured the risk score for every gene depending on their expression levels to determine the effects of the four IRGs on survival time and clinical prognosis in HCC patients, with the equation (Table 1):

In terms of the median measurement of risk, the HCC samples were divided into high- and low-risk categories. Then, we analyzed the impact of these risk groups on the overall survival proportion of HCC patients using KM curves. The high-risk group was found to have significantly worse OS (Figure 5A).

Figure 5:

Risk score analysis in TCGA. (A) Survival curve for low- and high-risk subgroups. (B) ROC curves for determining overall survival according to the risk score. (C) Heatmap showing gene expression, risk score distribution, and survival.

We conducted a time-dependent ROC analysis to evaluate the prognostic capability of the four IRGs. The area under the ROC curve (AUC) for the gene risk score signature was 0.721 (Figure 5B), showing moderate diagnostic capabilities. This suggests that the gene signature has the potential to be a reliable predictive tool. A heatmap was constructed to assess the levels of the four IRGs in the high and low-risk categories. Furthermore, the results of these four immunogenic biomarkers were analyzed in the ICGC-LIRI cohort. The results revealed that patients with an elevated risk measurement in the ICGC-LIRI cohort also showed poorer survival outcomes compared with those with a low-risk determination (Figure 6A–C). This finding suggests that the use of the four immune genes as biomarkers improved the sensitivity and specificity of the prognostic signature.

Figure 6:

Risk score assessment in ICGC. (A) Survival curve for low- and high-risk subgroups. (B) ROC curves for prediction of overall survival in terms of the risk score. (C) Heatmap showing gene expression, risk score distribution, and survival.

Construction of nomogram

According to the characteristics of the four immune-related genes, a nomogram was constructed to predict the survival outcomes of HCC patients (Figure 7). The estimated one-, three-, and five-year survival rates of HCC patients were determined using a vertical connection between the total point/each posterior axis. Furthermore, Cox regression analysis was used to evaluate the prognostic significance in relation to different clinical features. The tumor status, stage, T (tumor size), and risk score were found to be significantly linked with overall survival (p<0.01). Multivariate regression analysis showed that only the risk score was an independent predictor of survival outcome (p<0.01) (Table 2). This suggests that the risk score of the four IRGs provides valuable prognostic information, and acts as an independent predictor of patient outcomes. In summary, the nomogram incorporating the risk assessment with various clinical characteristics offers a comprehensive tool for forecasting the survival outcomes of HCC patients.

Figure 7:

Nomogram to determine one-, three-, and five-year overall survival using data from the TCGA cohort.

Validation of the predictive value

The KM-plotter tool was used for additional assessments. This tool determines the association between gene levels and OS. The log-rank test showed that the four IRGs were significantly associated with OS in HCC patients (p<0.01). This indicates that the levels of expression of these genes significantly influenced patient survival outcomes, supporting their prognostic potential as predictive biomarkers (Figure 8).

Figure 8:

Validation of the predictive value of IRGs using the Kaplan Meier-plotter in LIHC.

Immune status and GSEA of prognostic risk signature

We conducted PCA to compare the distribution patterns of the high and low-risk categories based on IRG expression and genome-level signaling. The IRG expression levels did not differ significantly from the immunological measurements based on the genome-level signaling features (Figure 9A). However, a significant difference was observed between the two groups, with lower scores seen in the group with minimal risk (Figure 9B).

Figure 9:

Immunological evaluation of the high- and low-risk groups using PCA and GSEA. (A) PCA of the high- and low-risk groups was performed using the whole protein-coding gene sets. (B) PCA of the high- and low-risk groups was performed using prognostic risk factors. (C and D) GO and KEGG enrichment analyses associated with risk levels.

Furthermore, GSEA was conducted on the high/low-risk signature datasets in TCGA-LIHC. GSEA indicated significant differences (NOM p<0.05) in the gene set enrichment from the MSigDB Collection. The most significantly enriched signaling bio-functions and network signals were identified following their normalized enrichment scores (NES). GSEA showed that genes were associated with modulation of the cell cycle transition between the G2 and M phases, spindle organization, complement activation, and humoral immune response (Figure 9C). The signaling pathways were primarily concentrated in base/nucleotide excision restoration, complement and coagulation cascades, and retinol metabolism (Figure 9D).

Correlation between the IRG prognostic signature and immune cell infiltration

Given the linkage between these IRGs and tumor immunity in HCC, data on immune cell infiltration were retrieved from the TIMER database. The analysis assessed the link between the IRG prognostic signature and immune cell infiltration in HCC. The results showed positive correlations between the risk measurement and infiltration of DCs, macrophages, neutrophils, and CD4+ T cells (Figure 10A–F). This suggests that the infiltration of these cell types is connected to HCC prognosis. The expression levels of the IRGs were associated with immune cell infiltration.

Figure 10:

Associations of risk factors and infiltration of six types of immune cells. Associations were determined using PCA. (A) DCs. (B) Macrophages. (C) Neutrophils. (D) CD4+ T cells. (E) CD8+ T cells. (F) B cells.