Sage Journals: Discover world-class research

Abstract

Objective:

This study aims to develop a prognostic model for HCC based on TME-related factors.

Introduction:

Hepatocellular carcinoma (HCC) is characterized by a poor prognosis, largely due to the complex and heterogeneous interactions between stromal and immune cells within the tumor microenvironment (TME).

Methods:

Genome and transcriptome data, as well as clinical information of HCC patients, were obtained from the Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO). The TME score was evaluated using the “ESTIMATE” R package. Differentially expressed genes (DEGs) associated with TME phenotype were analyzed using the LIMMA R-package. Survival outcomes were compared using Kaplan–Meier curves with log-rank test and Cox proportional hazards model. Protein-Protein Interaction (PPI) networks integrated with multivariate survival and LASSO analyses were utilized to identify TME-related hub genes for a risk score model. A nomogram predicting prognosis of HCC patients was developed through four independent cohorts.

Results:

The TME scores showed a negative correlation with tumor progression and survival in HCC patients. We identified 50 core genes with high connectivity in the PPI network, as along with 33 key DEGs associated with survival in HCC. Intersection analysis revealed six hub genes -CXCL8, CXCL1, CCR7, IL7R, MMP9, and CD69. The risk score based on these six TME-related hub genes was significantly associated with overall survival and clinicopathological characteristics of HCC patients. Furthermore, the nomogram demonstrated its ability to discriminate HCC patients from healthy individuals using peripheral blood mononuclear cells.

Conclusion:

We have developed a TME-related risk scoring model for HCC patients and identified six hub gene panel that serve as a potential biomarker for personalized prognosis of immunotherapy and non-invasive diagnostics of HCC.

Keywords

hepatocellular carcinoma tumor microenvironment differentially expressed genes risk score prognosis

Introduction

Hepatocellular carcinoma (HCC), the predominant histologic subtype of liver cancers, significantly contributes to global cancer-related morbidity and mortality.^1–3 It is also among the solid tumors with a poor prognosis. The incidence of HCC is steadily rising worldwide, leading to alarmingly high mortality rates. In 2022, approximately 865,000 new cases and 757,948 deaths caused by liver cancer were globally reported, resulting in a significant and challenging health burden.¹ Although chronic infection with hepatitis B virus (HBV) and hepatitis C virus (HCV) is accountable for the development of HCC, universally effective treatments for chronic infections with both HBV and HCV are contributing to the declines in the rates of viral-associated HCC. In contrast, the prevalence of metabolic risk factors is on the rise globally attributed to HCC, especially in developed countries. These factors include alcohol consumption, obesity, diabetes, and metabolic dysfunction-associated steatotic liver disease (MASLD), also known as non-alcoholic fatty liver disease (NAFLD).^1–5

HCC typically occurs in the presence of underlying cirrhosis and chronic liver inflammation.^2,3,6 A staggering 90% of HCC cases stem from cirrhosis, wherein hepatocytes undergo a relentless cycle of chronic necrosis and regeneration.⁶ Therefore, cirrhosis of the liver is the strongest risk factor for HCC.² HCV and HBV both cause chronic hepatitis and may lead to HCC through different mechanisms.^2,3 HCV infects hepatocytes, leading to activation of hepatic stellate cells and resulting in fibrosis. Its inflammatory response produces platelet-derived growth factors that further promote fibrosis.^2,6 Unlike HCV, HBV can integrate its genome with the host’s DNA, allowing its induced cases of hepatitis to progress to HCC even in the absence of fibrosis or cirrhosis.^2,7 Compared with viral-associated HCC, metabolic associated HCC presents a distinct molecular pathogenesis with genetic and nongenetic abnormalities,^2,3,6,7 highlighting the considerable heterogeneity in the development of HCC.

Environmental factors further exacerbate the risk of HCC. Alcoholic cirrhosis ranks as the second most common risk factor for HCC in the United States and Europe.^4,6 Studies have demonstrated a significant association between heavy alcohol consumption and an increased susceptibility to cirrhosis.⁸ Individuals with cirrhosis resulting from both hepatitis viral infection and alcohol face a significantly higher risk of HCC compared to those with alcohol-induced cirrhosis alone.⁹ The rates of HCV infection are 3–30 times higher in alcoholics than in the general population.¹⁰ HCV-infected patients who abuse alcohol exhibit more severe liver fibrosis and a heightened incidence of cirrhosis and HCC compared to non-alcohol-consuming counterparts.¹¹

Several recent studies have indicated that commonly used immune checkpoint inhibitors, such as nivolumab and pembrolizumab, have not shown statistically significant improvement in the primary endpoint of overall survival in phase III trials.^12,13 However, the combination of atezolizumab and bevacizumab (atezo-bev) has shown promising efficacy and safety as a first-line therapy for patients with unresectable hepatocellular carcinoma (uHCC) or Child-Turcotte-Pugh (CTP) class B cirrhosis.^14,15 Additionally, the benefits of tremelimumab in combination with durvalumab for a diverse global uHCC population were recently confirmed in a phase III clinical trial.^16,17 While tumor mutation burden (TMB) has been demonstrated to impact patient prognosis in many tumors,^18–20 it does not appear to affect patient survival in HCC due to its typically low TMB.²¹ The low remission rates of cancer immunotherapies for HCC patients can be attributed to the heterogeneity of the tumor microenvironment (TME) and immunosuppression.^22–24 Advancing safe and effective immunotherapies targeting TME is currently a key focus in cancer therapy.²⁵ TME is defined as the complex and diverse multicellular environment in which tumors develop. Interactions between malignant and non-transformed cells give rise to the TME, which plays an important role in malignant tumor progression, immune escape, and resistance to treatment.²⁵ Therefore, it is feasible to explore certain pivotal genes associate with the TME in HCC as prognostic indicators for HCC patients.

Given that stromal cells and immune cells are the two main non-tumor components of TME, they play a crucial role in cancer diagnosis and prognosis.^25,26 Therefore, utilizing the ESTIMATE algorithm proposed by Yoshihara to predict the levels of stromal and immune cell infiltration in TME is scientifically reliable.²⁷ Consequently, investigating the molecular basis of tumor immune microenvironment in HCC and establishing a risk score model could be beneficial for personalized immunotherapies.

The objective of his study was to establish a TME gene-related risk score for HCC prognosis using the ESTIMATE algorithm. We developed a nomogram model by integrating a six TME gene-based score with clinical features of HCC to provide precise prognostic insights. Validation through multiple datasets demonstrates that the hub gene panel (CXCL8, CXCL1, CCR7, IL7R, MMP9, and CD69) has discriminatory capacity for both prognosis and non-invasive diagnostics of HCC.

Materials and methods

Data sources

The RNA-seq gene expression data and corresponding clinical data of LIHC patients were obtained from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov).²⁸ A total of 423 RNA-seq data in TCGA-LIHC dataset, including 371 primary tumor tissues, 50 adjacent normal tissues and 2 metastatic tumor tissues, were utilized as the training cohort (Figure 1). The validation cohorts for this study were sought from various omics and diverse tissue types. The somatic mutation data of tumor DNA obtained from TCGA-LIHC were utilized as the validation cohort 1, providing validation from a distinct omics perspective. The validation cohort 2 comprised gene expression data of both HCC and non-cancerous tissues, measured by Affymetrix Human Genome U133 Plus 2.0 Array (GPL570), which were obtained from multiple datasets (GSE41804, GSE45267, GSE19665, GSE101685, and GSE121248) available on the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) repository.²⁹ The validation cohort 3 consisted of gene expression data of HCC measured by Affymetrix HT Human Genome U133A Array (GPL3921) and corresponding clinical prognosis information, which were obtained from GSE14520 dataset of GEO repository. The validation cohort 4 included gene expression data of peripheral blood mononuclear cell (PBMC) samples from both healthy individuals and patients with HCC, measured by GPL570 and obtained from GSE49515 dataset of GEO repository. The datasets for protein-protein interactions (PPIs) analysis were retrieved from the STRING database (https://string-db.org/).

Figure 1.

Flowchart of the study design.

This study used pre-defined datasets with established sample sizes, hence no adjustments were made to the criteria for including or excluding subjects. Detailed information about the samples is presented in Table S1, and the clinical characteristics of the HCC patients related to age, gender, fibrosis status, and risk factors was summarized in Table S2.

TME scores and their relationship with clinical characteristics

TME scores were estimated using the R package “ESTIMATE” from RNA-seq gene expression data for each sample in TCGA-LIHC dataset. The results yielded three scores: ImmuneScore, StromalScore, and EstimateScore, which represent the proportions of immune and stromal cells in the TME. A higher overall TME score indicates lower tumor purity. The comparison of TME scores was conducted between normal and primary tumor tissues. Additionally, further comparisons of TME scores were made between normal and primary tumor tissues stratified by different clinical stages and grades to explore their relationship with tumor progression.

To assess the correlation between TME scores and survival, patients in the top 25% and bottom 25% were categorized as the high- and low-score subgroups based on ImmuneScore and StromalScore, respectively. Subsequently, the Kaplan-Meier method was employed to compare the 5-year survival rates of these two groups and to generated survival curves.

The TME scores were compared between the group with history of alcohol consumption and the group without history of risk factors, the group with history of hepatitis and the group without ant known risk factors, respectively. Then survival analysis was performed between these groups.

DEG between high- and low-score groups and enrichment analysis

A total of 371 primary tumor samples were stratified into high- and low-score groups according to the median values of ImmuneScore and StromalScore, respectively. The “limma” package of R software was utilized to detect differentially expressed genes (DEGs) between high- and low-score groups. Genes with |log2FC|>1 and adjusted p-value < 0.05 were considered as DEGs. For ImmuneScore, genes with log2FC>1 and adjusted p-value < 0.05 were up-regulated genes in the high-score group, while genes with log2FC<-1 and adjusted p-value < 0.05 were down-regulated genes in high-score group. Similarly, distinct sets of up-regulated and down-regulated genes were also obtained for StromalScore. The DEGs associated with TME scores were determined by overlapping the ImmuneScore-related and StromalScore-related DEGs.

For TME scores-related DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed using “clusterProfiler” package of R software. Terms or pathways with an adjusted p-value < 0.05 were considered significantly enriched.

Protein-protein interaction (PPI) network

The PPI network was constructed using nodes (proteins) with a confidence score for their interactive relationships >0.7, and was subsequently visualized using Cytoscape software.³⁰ The degree of connectivity for each gene were assessed by counting the number of its neighboring nodes. The top 50 genes with the highest degree of connectivity were identified as the Core genes of the PPI network.

Survival-related genes in DEGs

Univariate survival analysis for each DEG was performed using the “survival” and “survminer” packages in R software. The 371 primary tumor samples were categorized into high- and low-expression groups based on the median TPM expression values of each gene, and Kaplan-Meier method with log-rank test was used to compare their survival rates. Gene expression levels were normalized using log2 (TPM+1), followed by univariate Cox proportional hazards regression analysis to calculate the p-value and HR (Hazard Ratio) at 95% confidence interval for each gene. A gene is deemed to an increased risk if its HR > 1, while a HR < 1 indicates a protective effect. A gene is deemed to be related to survival if both log-rank test and univariate Cox proportional hazards regression analysis yield p-values less than 0.05. The results of survival-related genes from the univariate survival analysis can be visualized in a forest plot.

Hub genes and multi-variates cox analysis

By intersecting the Core genes in the PPI network with the survival-related genes in the DEGs, a set of pivotal genes known as hub genes were identified. These hub genes play a crucial role within the PPI network and demonstrate associations with both TME scores and survival outcomes. To identify independent prognostic factors, multivariate Cox proportional hazards regression analysis was performed on the hub genes, incorporating other clinical characteristics including gender, clinical stage and grade information, using the “survival” and “survminer” packages in R software.

Risk score model and nomogram

To eliminate the possible gene collinearity and obtain the optimal combination of genes, Lasso Cox analysis on hub genes was conducted using the “glmnet” and “survival” packages in R software.³¹ The λ value corresponding to the minimum partial likelihood deviance, determining through 10-fold cross validations, was selected as the optimal λ for our study. A risk score for survival prediction was calculated based on TCGA-LIHC gene expression data in the training cohort using the formula: Risk score = $\sum β_{i} E_{i}$ , where $i \in$ hub genes, $E_{i}$ is the log2(TPM+1) expression value of the gene $i,$ and $β_{i}$ denotes the coefficients of gene $i$ in Lasso Cox regression.

TCGA-LIHC patients were stratified into high- and low-risk groups based on the median risk score, and their survival was analyzed using Kaplan-Meier method. Multivariate Cox proportional hazards regression analysis was performed to assess whether the risk score served as an independent prognostic factor. Additionally, correlation analyses were conducted between patients’ risk score and their TME-scores, as well as between risk score and survival time for deceased patients. Furthermore, comparisons of risk scores were made between normal and primary tumor tissues categorized by different clinical stages and grades.

Nomogram analysis was performed in the TCGA-LIHC dataset using “rms” package in R. The total points and sum of points for every factor could accurately predict the 1-, 3-, and 5-year survival rates of patients.

Verifications and validations

To enhance the credibility of the risk score model of six hub genes, multiple datasets were employed for validation. The TCGA-LIHC somatic mutation data from tumor DNA (validation cohort 1) was used to stratify TCGA-LIHC patients into high- and low-risk groups based on the median risk score, and then compared TMB of each patient between these two groups. Fisher exact test was further performed to identify differentially mutated genes with adjusted p-value < 0.05 between low- and high-risk groups. Using the gene expression data from GSE41804, GSE45267, GSE19665, GSE101685, and GSE121248 datasets (validation cohort 2), the expressions of each hub gene were compared between normal and HCC tissues. The gene expression data and corresponding clinical prognosis information from GSE14520 dataset (validation cohort 3) were utilized for survival analysis using univariate Cox regression model for each hub gene and the risk score. The potential of hub genes for the non-invasive testing were assessed using principal component analysis (PCA) of the gene expression data from PBMC samples in GSE49515 dataset, including both healthy individuals and patients with HCC (validation cohort 4).

Statistical analysis

The comparison of values between two groups was analyzed using the Wilcoxon rank sum test, and the relationship between pairs of values was evaluated through Spearman correlation analysis. Statistical signiﬁcance was defined as p < 0.05, denoted as *, p < 0.01 as **, and p < 0.001 as ***.

Results

TME scores associated with tumor progression and prognosis of HCC patients

The patients from TCGA-LIHC dataset (Table S1) were initially divided into high- and low-score groups based on their top 25% and bottom 25% of TME scores, respectively. Significantly lower TME scores were observed in HCC tissues (n = 371) compared to the normal controls (n = 50; p < 0.001; Figure 2(a)), even in their paired tissues (n = 50; p < 0.001; Figure 2(b)). The TME scores exhibited a gradual decrease in correlation with the progression of tumor stages (Figure 2(c)) and pathological grades (Figure 2(d)), except for stage IV (n = 5) and grade IV (n = 12). These data suggest an inverse relationship between the TME scores and tumor progression in HCC patients (p < 0.001). Furthermore, the 5-year survival of the patients was significantly associated with improved survival time based on Stromal scores (p < 0.05), Immune scores (p < 0.01) and ESTIMATE scores (p < 0.05; Figure 2(e)), indicating that higher TME scores may serve as a predictor for better prognosis in HCC.

Figure 2.

The clinical relevance of TME scores in HCC patients. (a–b) Comparison of TME scores between HCC tumors and non-cancerous tissues. (a) HCC tumors (n = 371) versus noncancerous tissues (Normal, n = 50). ***p < 0.001, by Wilcoxon rank sum test. (b) Paired tissue (n = 50) (c–d) Variation in TME scores among patients at different tumor stages (c) and with different pathological grades (d), compared to normal controls. *p < 0.05, **p < 0.01, p < 0.001, by Wilcoxon rank sum test. (e) The 371 patients from the TCGA-LIHC dataset were divided into high- and low-score groups based on their top 25% and bottom 25% of TME scores, respectively. Kaplan-Meier curves were used to analyze the 5-year survival of the patient between high- and low-score groups.

Impacts of alcohol consumption and hepatitis on TME in HCC

Alcohol consumption and hepatitis are two of the most common clinical risk factors for HCC, leading to chronic inflammation. In the TCGA-LIHC dataset, 118 patients had alcohol consumption, 114 had hepatitis, and 92 had no history of primary risk factors (Table S2). Patients with alcohol consumption showed significantly higher Immune scores (p<0.01; Figure 3(a)), while those with hepatitis showed slightly higher but not statistically significant TME scores (Figure 3(b); Table S3). HCC patients with alcohol consumption or hepatitis had better 5-year survival than those without history of primary risk factors (p<0.05; Figure 3(c)). When comparing the significant TME scores of patients with combined factors of cirrhosis, age, and gender, HCC patients with established cirrhosis only demonstrated significantly higher Immune score (Tables S3 and S4), confirming the immune impact of TME on HCC.

Figure 3.

Impacts of alcohol consumption and hepatitis on TME scores in HCC. HCC patients from the TCGA-LIHC cohort divided into three groups: those with alcohol consumption (n = 117), those with hepatitis (n = 114) and those without history of primary risk factors (n = 91). (a–b) The TME scores of HCC patients without history of primary risk factor were compared to those with alcohol consumption (a), or hepatitis (b), using the Wilcoxon rank sum test for p value analysis. (c) The 5-year survival rates of patients with alcohol consumption or hepatitis were compared to those without a history of primary risk factors.

TME-related gene in HCC

To study TME-related genes expressed in HCC, we divided the TCGA-LIHC tumor samples into high-score (n = 185) low-score (n = 186) groups based on their median TME scores. The heatmaps in Figure 4(a) show gene expression patterns related to TME score in HCC. In the high-score group, we identified 1338 up-regulated and 84 down-regulated genes associated with stromal scores (Figure S1A), as well as 1041 up-regulated and 81 down-regulated gene linked to immune scores (Figure S1B). A total of 830 TME-related differentially expressed genes (DEGs) were selected through overlapping between stromal and immune scores, consisting of 802 up-regulated and 28 down-regulated DEGs (Figure 4(b)). These TME-related DEGs were primarily enriched in terms such as leukocyte activation and adhesion (Figure 4(c)), with top functional categories including T cell activation, leukocyte cell-cell adhesion, leukocyte proliferation, regulation of cell-cell adhesion, and positive regulation of cell adhesion (Figure 4(d)). Additionally, KEGG pathway analysis showed that these DEGs were mainly enriched in pathways such as cytokine-cytokine receptor interaction, chemokine signaling pathway, and P3K-Akt signaling pathway (Figure 4(d)).

Figure 4.

TME score-related DEGs in HCC. A total of 371 TCGA-LIHC patients were categorized into high- and low-score groups based on the median value of Immune score and Stromal score, respectively. These DEGs were considered to be relevant to the TME. (a) The heatmaps showing association between the expression of DEGs and Immune score (Up panel) as well as stromal score (Down panel). (b) Among the overlapped DEGs related to Immune score- or Stromal score in HCC, 802 were up-regulated (Left), while 28 were down-regulated DEGs (Right). (c) The top five categories by GO enrichment of the TME-related DEGs including leukocyte-cell adhesion, regulation of leukocyte-cell adhesion, regulation of lymphocyte activation, regulation of T-cell activation, T-cell activation, and T-cell activation. Detailed information about these genes is provided. (d) GO (Left) and KEGG enrichment (Right) analysis of the TME-related DEGs, respectively.

Identification of TME-related hub genes in HCC

The PPI network of the TME-related DEGs consists of 544 nodes and 2587 edges, representing confident interactions (Figure 5(a)). The top 50 DEGs with largest degrees were identified as the core genes in the PPI network (Figure 5(b)). Subsequently, univariate cox proportional hazards regression and log-rank test were performed on 830 TME-related DEGs, leading to the identification of 33 DEGs associated with HCC prognosis (Figure 5(c)). Furthermore, through interaction analysis between the core genes and the prognosis-related DEGs in HCC, a panel of six genes, including CXCL8, CXCL1, CCR7, IL7R, MMP9, and CD69, was identified (Figure 5(d)). Each of these six genes is significantly related with the survival of HCC (Figure S1C) and thus considered as hub genes. In a training cohort from TCGA-LIHC (Figure 1), multivariate cox regression analyses showed that inflammatory chemokine CXCL8 independently contributes to HCC prognosis (Figure 5(e)).

Figure 5.

Identification of hub genes related with TME scores. (a) The PPI network of TME-related DEGs in HCC. (b) The top 50 of TME-related DEGs with the largest degrees, serving as the Core genes in the PPI network. (c) Forest plot of univariate Cox regression analysis of survival-related genes among TME-related DEGs, showing p values and hazard ratio (HR) with 95% confidence interval (CI) for overall survival of DEGs. (d) Six Hub genes identified through the overlap of the 50 Core genes and 33 survival-related DEGs. (e) Forest plot of multivariate Cox regression analysis for each individual Hub gene.

Construction of TME gene-based risk score model for HCC prognosis

A TME gene-based risk score model for HCC prognosis was developed with LASSO Cox regression of the six hub genes. The model with the optimal parameter value (λ = 0.007) was selected based on the minimum partial likelihood deviance through 10-fold cross validations (Figure 6(a)). All six hub genes based on the optimal λ were retained in the final risk score model. For each sample, the risk score can be calculated as $\begin{array}{l} E_{C X C L 8} \times 0.279 + E_{C X C L 1} \times (- 0.06) \\ + E_{I L T R} \times (- 0.378) + E_{M M P 9} \times 0.132 \\ + E_{C C R 7} \times (- 0.072) + E_{C D 69} \times (- 0.027), \end{array}$ where $E_{i}$ represents the log2(TPM+1) expression value of the gene $i,$ if gene expression values are represented differently, adjustments should be made the coefficients accordingly.

Figure 6.

Hub gene-based risk score model for HCC prognosis. (a) Construction of LASSO Cox regression with six hub genes. Left: the curves show the change track of each independent variable coefficient of the hub genes (Y-axis) in relation to the parameter lambda (λ) (X-axis). Each of the hub genes are shown with different colors. Right: the partial likelihood deviance respect to parameter λ. (b) Comparison of risk scores between the HCC patients (n = 371) and normal controls (n = 50). (c) Correlation analysis between risk score and TME scores of the HCC patients (n = 371). Left: with Stromal score. Right: with Immune score. (d) Multivariate Cox regression was performed to assess the independent prognostic value of the risk score. (e) Based on the median of the risk score, 371 TCGA-LIHC patients were assigned into high- and low-risk score groups for survival analysis using Kaplan-Meier curves. (f) Spearman correlation analysis of the relationship between the risk score and survival time in LIHC-HCC patients with the deceased status. (g) A nomogram model integrating the six TME gene-based risk score with gender, alcohol consumption and hepatitis was developed to predict the 1-, 3-, and 5-year survival rates of HCC patients from TCGA-LIHC cohort. (h) The receiver operating characteristic (ROC) curve showed higher area under the curves (AUCs), specifically 0.768 for 1-year, 0.728 for 2-year, and 0.751 for a 3-year survival rate for based on nomogram six hub genes in TCGA training cohort.

The TCGA-LIHC cohort was utilized to assess the clinical relevance of the risk score in HCC. In HCC tissues (n = 371), a higher risk score was observed compared to non-cancerous tissues (n = 50; p < 0.01; Figure 6(b)). Spearman correlation analysis revealed a negative correlation between the risk score and TME scores of the HCC patients (Figure 6(c)). Additionally, multivariate cox regression analysis of TCGA-LIHC cohort indicated that the risk score was an independent prognostic factor for survival of HCC patients (p < 0.001; Figure 6(d)). The 371 HCC tumor samples were thus categorized into high-risk or low-risk groups, with survival analysis showing inferior outcomes for patients in the high-risk group compared to those in the low-risk group (p < 0.001; Figure 6(e)). Among deceased HCC patients, a negative correlation was observed between risk score and overall survival time (p = 0.01; Figure 6(f)). These findings suggest that the TME gene-based risk score assessment is informative for predicting survival outcomes in HCC.

A nomogram was developed in the TCGA-LIHC cohort to establish a practical method for clinical predicting the prognosis of HCC patients. This nomogram integrates gender, alcohol consumption, hepatitis, and six TME gene-based risk score for 1-, 3-, and 5-year survival prediction. Figure 6(g) illustrates the accurate conversion between total points and the probability of clinical outcomes. The effectiveness of the nomogram model on HCC prognosis was validated by evaluating and comparing the accuracy of 1-, 2-, and 3-year survival using receiver operating characteristic (ROC) analysis, which yielded values of 0.768 for 1-year, 0.728 for 2-year, and 0.751 for 3-year (Figure 6(h)).

Multi-cohorts verifications and validations

To further investigate the impact of genetic factors on the risk score, TCGA-LIHC patients with tumor DNA somatic mutation data (n = 362, validation cohort 1) were stratified into high- and low-risk groups. It was observed that the low-risk group exhibited significantly lower TMBs compared to the high-risk group (p < 0.05; Figure 7(a)). Additionally, the TP53 mutation rate was found to be significantly higher in the high-risk group than in the low-risk group, with TP53 being identified as the only gene showing significant differential mutation between both groups (p < 0.001; Figure 7(b), Figure S2A). Furthermore, a combined cohort consisting 286 samples (validation cohort 2), included 172 HCC tumor tissues and 114 normal controls, was utilized for validating the expression of 6 TME hub genes. The results demonstrated that mRNA levels of five out of six hub gene was significantly downregulated, while MMP9 showed upregulation (Figure 7(c)), thus confirming their differentially expressed in HCC.

Figure 7.

Verifications and validations of the TME score-related gene panel using multi-datasets. (a–b) In validation cohort 1 (n = 362) from TCGA-LIHC, tumor DNA somatic mutation data were used to compare high- and low-risk groups to assess genetic factors contributing to the risk score. Tumor mutation burdens were analyzed by the Wilcoxon rank sum test (*p < 0.05) (a) and TP53 status was assessed by Fisher exact test (***p < 0.001) (B). (C) Validation cohort 2 (n = 286), formed by combing GES 41804, GSE45267, GSE19665, GSE101685, and GES121248 datasets, included 172 HCC tumor tissues and 114 normal controls. The differential expression of the six hub genes between normal and HCC tissue samples was analyzed using the Wilcoxon rank sum test. (d) GES14520, referred to as validation cohort 3 (n = 225), was utilized for univariate Cox regression to access the clinical relevance of the TME score-related gene panel risk score. The forest plot illustrates the association of the risk score with overall survival and recurrence in HCC patients. (e) Peripheral blood mononuclear cell (PBMC) samples from GSE49515 (10 HCC patients and 10 healthy individuals) were recruited as the validation cohort 4. Discrimination between HCC patients and 10 healthy individuals was achieved by using DEGs, including a set of over 23,000 global DEGs (Left panel), or the 6 TME score-related gene panel (Right panel).

Next, the forest plot of univariate cox regression on the GSE14520 dataset (validation cohort 3, n = 225) reveals a significant association between the hub gene-based risk score and both survival and recurrence of HCC. However, it is important to note that a single hub gene may not possess significant prognostic value (Figure 7(d)). Furthermore, gene expression data from PBMC samples was included for healthy individuals (n = 10) and patients with HCC (n = 10) (validation cohort 4). The principal component analysis demonstrates that healthy individuals and HCC patients can be distinguished with either more than 23,000 genes or the six-hub gene panel (Figure 7(e)). These findings validate the prognostic significance of the six TME-gene panel and its potential for noninvasive testing.

Discussion

Numerous studies have confirmed an essential role of TME in tumor growth, invasion, and metastasis, as well as its association with prognosis and therapy response.^{25,26,32–34} The survival analysis in this study further confirmed that a high TME score is associated with a favorable survival outcome in HCC. A high TME score often indicates a more favorable microenvironment for anti-tumor responses, characterized by higher immune and stromal cell infiltration, higher immunogenicity, and lower tumor purity.³⁵ Conversely, a low TME score may indicate a less responsive microenvironment with limited immune and stromal cell infiltration, potentially promoting tumor progression. Significant differences in TME scores were also observed between normal tissues and HCC of distinct stages and grades. It is noteworthy that patients with alcohol consumption or hepatitis tends to have higher Immune scores and better survival compared to those without history of primary risk factors. Further investigation into their genetic backgrounds revealed that these patients with a history of alcohol consumption or hepatitis exhibits higher somatic mutation burdens (Figure S2B). Some studies^36,37 suggest that high somatic mutation burden can result in the expression of neoantigens on cancer cells, potentially rending the tumor more visible to the immune system, increasing the likelihood of an anti-tumor immune response.

Although TME scores have proved valuable information, the acquisition of bulk gene expression data for calculating TME scores by the ESTIMATE algorithm³⁸ can be costly. To address this limitation, we have identified six hub genes (CXCL8, CXCL1, CCR7, IL7R, MMP9, and CD69) that are associated with both TME and the prognosis of HCC. This was achieved through a combination of correlation analysis and survival analysis for the TME-related DEGs. Each of these hub genes demonstrated significant differences in expression between HCC and non-cancerous tissues. Enrichment analysis revealed that these genes primarily functions in tumor immunology, which aligns with their relevance to the TME. Among them, CXCL8 and CXCL1 are two significant chemokines that facilitate tumor cell proliferation, migration, invasion, and contribute to the formation of the TME.³⁹ CXCL8, also known as IL-8, is a key mediator of the inflammatory response and plays a vital role in the invasion and metastasis of HCC. It is closely associated with poor prognosis of HCC patients.⁴⁰ Similarly, CXCL1, a commonly elevated inflammatory chemokine in HCC, directly enhances tumorigenesis within the hepatic inflammatory microenvironment by acting on tumor cells and indirectly through promoting macrophages recruitment.⁴¹ CCR7 serves as the chemokine receptor of CCL19 and CCL21, whose interaction regulates tumor cell migration and proliferation, or controls the recruitment of tumor-related microglia/macrophage.⁴² IL7R, as the receptor of IL-7, may be involved in regulating lymphocyte infiltration within the TME, thereby influencing the progression of tumor.⁴³ The expression of MMP9 is associated with regulatory T-cells infiltration in the TME of HCC.⁴⁴ In the context of anti-tumor immunity, CD69 negatively regulates the effector function of intratumoral T cells and importantly controls the exhaustion of CD8 T cells.⁴⁵ These hub genes can serve as reliable biomarkers, and their differential expression levels in HCC can provide valuable information about the TME and its impact on disease progression.

We have developed a prognosis model using LASSO Cox regression, utilizing only six hub genes to calculate a risk score for each patient. The risk score demonstrated a negative correlation with both Immune score and Stromal score, indicating its relevance to the TME. Survival analysis revealed that risk score is an independent prognostic factor in HCC. In both the training cohort and validation cohort, a high-risk score is associated with unfavorable survival outcomes in HCC. This risk score could serve as a robust prognostic biomarker for predicting survival and recurrence risk in HCC patients, guiding healthcare professionals in precise treatment and rehabilitation strategies. Further comparison of genetic backgrounds revealed that patients in the high-risk group exhibits higher somatic mutation burdens than their counterparts in the low-risk group. Additionally, the mutation rate of TP53 was significantly higher in the high-risk group compared to the low-risk group. As a well-known tumor suppressor gene, both TP53 missense and nonsense mutations are related with elevated TMB and neoantigen levels, contributing to deficiencies in DNA damage repair.⁴⁶ The hub genes and TP53 were involved multiple GO categories and tumor immune pathways (Figure S2C). As a novel prognosis biomarker based on TME of HCC, the risk score also demonstrates significant non-invasive diagnostic potential, effectively distinguishing between normal individuals and HCC patients.

The primary contribution of this study is the introduction of a risk score model based on six hub genes, which is capable of generating a personalized risk score for each patient. The risk score serves as an independent prognostic factor for HCC and holds significant potential for broad applications in clinical practice. However, the study has certain limitations. First, this study relied on the analysis of multiple publicly available datasets for validation without additional sample size calculation and did not employ our own samples for validation. Second, the current risk model was based on TPM values of RNA-seq; therefore, the coefficients may require further adjustment at protein levels in the future, particularly about the design of this six-gene panel.

Conclusions

In this study, a nomogram model integrating a score based on six TME genes along with the clinical features of HCC was developed. This six TME gene panel (CXCL8, CXCL1, CCR7, IL7R, MMP9, and CD69) might possess discriminatory capabilities for both the prognosis and non-invasive diagnostics of HCC.

Supplemental Material

sj-docx-2-iji-10.1177_03946320251333975 – Supplemental material for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma

Supplemental material, sj-docx-2-iji-10.1177_03946320251333975 for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma by Xinyi Li, Zifan Qin, Haozhi Chen, Daichuan Chen, Nafisa Alimu, Duoduo Li, Xiyu Cheng, Qiong Yan, Lishu Zhang, Xingwei Liu, Zitong Zhou, Jiayi Zhu, Hangqi Ma, Xinyue Pei, Hanli Xu and Jiaqiang Huang in International Journal of Immunopathology and Pharmacology

Supplemental Material

sj-docx-3-iji-10.1177_03946320251333975 – Supplemental material for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma

Supplemental material, sj-docx-3-iji-10.1177_03946320251333975 for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma by Xinyi Li, Zifan Qin, Haozhi Chen, Daichuan Chen, Nafisa Alimu, Duoduo Li, Xiyu Cheng, Qiong Yan, Lishu Zhang, Xingwei Liu, Zitong Zhou, Jiayi Zhu, Hangqi Ma, Xinyue Pei, Hanli Xu and Jiaqiang Huang in International Journal of Immunopathology and Pharmacology

Supplemental Material

sj-docx-4-iji-10.1177_03946320251333975 – Supplemental material for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma

Supplemental material, sj-docx-4-iji-10.1177_03946320251333975 for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma by Xinyi Li, Zifan Qin, Haozhi Chen, Daichuan Chen, Nafisa Alimu, Duoduo Li, Xiyu Cheng, Qiong Yan, Lishu Zhang, Xingwei Liu, Zitong Zhou, Jiayi Zhu, Hangqi Ma, Xinyue Pei, Hanli Xu and Jiaqiang Huang in International Journal of Immunopathology and Pharmacology

Supplemental Material

sj-jpg-5-iji-10.1177_03946320251333975 – Supplemental material for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma

Supplemental material, sj-jpg-5-iji-10.1177_03946320251333975 for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma by Xinyi Li, Zifan Qin, Haozhi Chen, Daichuan Chen, Nafisa Alimu, Duoduo Li, Xiyu Cheng, Qiong Yan, Lishu Zhang, Xingwei Liu, Zitong Zhou, Jiayi Zhu, Hangqi Ma, Xinyue Pei, Hanli Xu and Jiaqiang Huang in International Journal of Immunopathology and Pharmacology

Supplemental Material

sj-jpg-6-iji-10.1177_03946320251333975 – Supplemental material for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma

Supplemental material, sj-jpg-6-iji-10.1177_03946320251333975 for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma by Xinyi Li, Zifan Qin, Haozhi Chen, Daichuan Chen, Nafisa Alimu, Duoduo Li, Xiyu Cheng, Qiong Yan, Lishu Zhang, Xingwei Liu, Zitong Zhou, Jiayi Zhu, Hangqi Ma, Xinyue Pei, Hanli Xu and Jiaqiang Huang in International Journal of Immunopathology and Pharmacology

Supplemental Material

sj-xlsx-1-iji-10.1177_03946320251333975 – Supplemental material for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma

Supplemental material, sj-xlsx-1-iji-10.1177_03946320251333975 for Construction of a tumor immune microenvironment-related risk scoring model for prognosis of hepatocellular carcinoma by Xinyi Li, Zifan Qin, Haozhi Chen, Daichuan Chen, Nafisa Alimu, Duoduo Li, Xiyu Cheng, Qiong Yan, Lishu Zhang, Xingwei Liu, Zitong Zhou, Jiayi Zhu, Hangqi Ma, Xinyue Pei, Hanli Xu and Jiaqiang Huang in International Journal of Immunopathology and Pharmacology

Footnotes

Acknowledgements

We would like to express our heartfelt gratitude to Mr. Enkai Wang, Mr. Qihang Peng, and all our laboratory members for their generous assistance with the manuscript revision.

Authors’ contribution

(1) H. X, J. H, X. L: Conception or design; X. L, Z. Q, H. C, H. X, J. H: Methodology: X. C, Q. Y, L. Z; Resources: D. C, X. L, Z. Z; J. Z, H. M, X. P, N.A, D.L; Writing—original draft preparation: X. L, Z. Q, H. C.

(2) Writing—review and editing: X. L, H, C, H. X, J. H; Supervision: H. X, J. H.

(3) Final approval of the manuscript: All authors.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was funded by The JIAODA-NORSON Innovation Research and Development Program (KSL18030530); Beijing Natural Science Foundation (7242091); Beijing Jiaotong University undergraduate innovation and entrepreneurship training project (No. 20231000411012, 2023100041500, 20231000411295).

Ethics approval

Not applicable.

Informed consent

Not applicable.

Trial registration

Not applicable.

ORCID iD

Jiaqiang Huang

Data availability

Not applicable.

Code availability

Not applicable.

Supplemental material

Supplemental material for this article is available online.

References

Bray

Laversanne

Sung

, et al. (2024) Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: A Cancer Journal for Clinicians 74: 229–263.

Llovet

Kelley

Villanueva

, et al. (2021) Hepatocellular carcinoma. Nature Reviews Disease Primers 7: 6.

Toh

Wong

EYT

Wong

, et al. (2023) Global epidemiology and genetics of hepatocellular carcinoma. Gastroenterology 164: 766–782.

Kim

(2024) Changing etiology and epidemiology of hepatocellular carcinoma: Asia and worldwide. Journal of Liver Cancer 24: 62–70.

Hanbury

Jani

Abdallah

, et al. (2023) Geographic and temporal trends in aetiology, incidence and mortality from hepatocellular carcinoma in European Union 15+ countries. Liver Cancer International 4: 109–120.

Kim

Viatour

(2020) Hepatocellular carcinoma: Old friends and new tricks. Experimental & Molecular Medicine 52: 1898–1907.

Armandi

Rosso

Caviglia

, et al. (2025) An updated overview on hepatocellular carcinoma in patients with metabolic dysfunction-associated steatotic liver disease: Trends, pathophysiology and risk-based surveillance. Metabolism 162: 156080.

Huang

Mathurin

Cortez-Pinto

, et al. (2023) Global epidemiology of alcohol-associated cirrhosis and HCC: Trends, projections and risk factors. Nature Reviews Gastroenterology & Hepatology 20: 37–49.

Abassa

Xiao

, et al. (2022) Effect of alcohol on clinical complications of hepatitis virus-induced liver cirrhosis: A consecutive ten-year study. BMC Gastroenterology 22: 130.

10.

Kouroumalis

Voumvouraki

(2022) Hepatitis C virus: A critical approach to who really needs treatment. World Journal of Hepatology 14: 1–44.

11.

Pearson

Kim

Berry

, et al. (2021) Associations between alcohol use and liver-related outcomes in a large national cohort of patients with cirrhosis. Hepatology Communications 5: 2080–2095.

12.

Chang Lee

Tebbutt

(2019) Systemic treatment of advanced hepatocellular cancer: New hope on the horizon. Expert Review of Anticancer Therapy 19: 343–353.

13.

Liu

Lin

Zhang

, et al. (2019) Molecular targeted and immune checkpoint therapy for advanced hepatocellular carcinoma. Journal of Experimental & Clinical Cancer Research 38: 447.

14.

Kulkarni

Tevethia

Kumar

, et al. (2023) Effectiveness and safety of atezolizumab-bevacizumab in patients with unresectable hepatocellular carcinoma: A systematic review and meta-analysis. EClinicalMedicine 63: 102179.

15.

Song

Yeom

Jung

, et al. (2024) Risk of bleeding in hepatocellular carcinoma patients treated with atezolizumab/bevacizumab: A systematic review and meta-analysis. Liver Cancer 13: 590–600.

16.

Sangro

Chan

Kelley

, et al. (2024) Four-year overall survival update from the phase III HIMALAYA study of tremelimumab plus durvalumab in unresectable hepatocellular carcinoma. Annals of Oncology 35: 448–457.

17.

Abou-Alfa

Lau

Kudo

, et al. (2022) Tremelimumab plus durvalumab in unresectable hepatocellular carcinoma. NEJM Evidence 1: EVIDoa2100070.

18.

Samstein

Lee

Shoushtari

, et al. (2019) Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nature Genetics 51: 202–206.

19.

Wang

Liu

Yang

, et al. (2022) A novel tumor mutational burden-based risk model predicts prognosis and correlates with immune infiltration in ovarian cancer. Frontiers in Immunology 13: 943389.

20.

Yarchoan

Hopkins

Jaffee

(2017) Tumor mutational burden and response rate to PD-1 inhibition. The New England Journal of Medicine 377: 2500–2501.

21.

Mauriello

Zeuli

Cavalluzzo

, et al. (2019) High somatic mutation and neoantigen burden do not correlate with decreased progression-free survival in HCC patients not undergoing immunotherapy. Cancers 11: 1824.

22.

Filliol

Saito

Nair

, et al. (2022) Opposing roles of hepatic stellate cell subpopulations in hepatocarcinogenesis. Nature 610: 356–365.

23.

Yang

Quan

, et al. (2022) A single-cell atlas of the multicellular ecosystem of primary and metastatic hepatocellular carcinoma. Nature Communications 13: 4594.

24.

Zhu

Finn

Edeline

, et al. (2018) Pembrolizumab in patients with advanced hepatocellular carcinoma previously treated with sorafenib (KEYNOTE-224): A non-randomised, open-label phase 2 trial. The Lancet. Oncology 19: 940–952.

25.

Bejarano

Jordao

MJC

Joyce

(2021) Therapeutic targeting of the tumor microenvironment. Cancer Discovery 11: 933–959.

26.

Mun

Leem

Lee

, et al. (2022) Dual relationship between stromal cells and immune cells in the tumor microenvironment. Frontiers in Immunology 13: 864739.

27.

Yoshihara

Shahmoradgoli

Martinez

, et al. (2013) Inferring tumour purity and stromal and immune cell admixture from expression data. Nature Communications 4: 2612.

28.

Cancer Genome Atlas Research Network; Weinstein

Collisson

, et al. (2013) The Cancer Genome Atlas Pan-Cancer analysis project. Nature Genetics 45: 1113–1120.

29.

Edgar

Domrachev

Lash

(2002) Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Research 30: 207–210.

30.

Shannon

Markiel

Ozier

, et al. (2003) Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Research 13: 2498–2504.

31.

Friedman

Hastie

Tibshirani

(2010) Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software 33: 1–22.

32.

Quail

Joyce

(2013) Microenvironmental regulation of tumor progression and metastasis. Nature Medicine 19: 1423–1437.

33.

Tiwari

Trivedi

Lin

(2022) Tumor microenvironment: Barrier or opportunity towards effective cancer therapy. Journal of Biomedical Science 29: 83.

34.

Hinshaw

Shevde

(2019) The tumor microenvironment innately modulates cancer progression. Cancer Research 79: 4557–4566.

35.

Zhang

, et al. (2021) A risk model developed based on tumor microenvironment predicts overall survival and associates with tumor immunity of patients with lung adenocarcinoma. Oncogene 40: 4413–4424.

36.

Jardim

Goodman

de Melo Gagliato

, et al. (2021) The challenges of tumor mutational burden as an immunotherapy biomarker. Cancer Cell 39: 154–173.

37.

Liu

Tan

, et al. (2021) High-affinity neoantigens correlate with better prognosis and trigger potent antihepatocellular carcinoma (HCC) activity by activating CD39⁺CD8⁺ T cells. Gut 70: 1965–1977.

38.

Yang

Liu

Nan

, et al. (2020) Comprehensive analysis of prognostic immune-related genes in the tumor microenvironment of cutaneous melanoma. Journal of Cellular Physiology 235: 1025–1035.

39.

Lin

Zhou

, et al. (2023) New genetic and epigenetic insights into the chemokine system: The latest discoveries aiding progression toward precision medicine. Cellular & Molecular Immunology 20: 739–776.

40.

Sun

Wang

Sun

, et al. (2019) Interleukin-8 promotes integrin β3 upregulation and cell invasion through PI3K/Akt pathway in hepatocellular carcinoma. Journal of Experimental & Clinical Cancer Research 38: 449.

41.

Zhao

Wei

Zhou

, et al. (2023) Blocking the CXCL1-CXCR2 axis enhances the effects of doxorubicin in HCC by remodelling the tumour microenvironment via the NF-κB/IL-1β/CXCL1 signalling pathway. Cell Death Discovery 9: 120.

42.

Geraldo

Garcia

, et al. (2023) CCL21-CCR7 signaling promotes microglia/macrophage recruitment and chemotherapy resistance in glioblastoma. Cellular and Molecular Life Sciences 80: 179.

43.

Wang

Chang

Wang

, et al. (2022) IL7R is correlated with immune cell infiltration in the tumor microenvironment of lung adenocarcinoma. Frontiers in Pharmacology 13: 857289.

44.

Zhang

(2022) Regulatory T-cells-related signature for identifying a prognostic subtype of hepatocellular carcinoma with an exhausted tumor microenvironment. Frontiers in Immunology 13: 975762.

45.

Koyama-Nasu

Wang

Hasegawa

, et al. (2022) The cellular and molecular basis of CD69 function in anti-tumor immunity. International Immunology 34: 555–561.

46.

Sun

Liu

Zhou

, et al. (2020) Specific TP53 subtype as biomarker for immune checkpoint inhibitors in lung adenocarcinoma. EBioMedicine 60: 102990.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.02 MB

1.55 MB

0.18 MB