Identification of Prognostic Genes in the Tumor Microenvironment of Hepatocellular Carcinoma.

Background: Hepatocellular carcinoma (HCC) is one of the most common malignant tumors in the world. The efficacy of immunotherapy usually depends on the interaction of immunomodulation in the tumor microenvironment (TME). This study aimed to explore the potential stromal-immune score-based prognostic genes related to immunotherapy in HCC through bioinformatics analysis.
Methods: ESTIMATE algorithm was applied to calculate the immune/stromal/Estimate scores and tumor purity of HCC using the Cancer Genome Atlas (TCGA) transcriptome data. Functional enrichment analysis of differentially expressed genes (DEGs) was analyzed by the Database for Annotation, Visualization, and Integrated Discovery database (DAVID). Univariate and multivariate Cox regression analysis and least absolute shrinkage and selection operator (LASSO) regression analysis were performed for prognostic gene screening. The expression and prognostic value of these genes were further verified by KM-plotter database and the Human Protein Atlas (HPA) database. The correlation of the selected genes and the immune cell infiltration were analyzed by single sample gene set enrichment analysis (ssGSEA) algorithm and Tumor Immune Estimation Resource (TIMER).
Results: Data analysis revealed that higher immune/stromal/Estimate scores were significantly associated with better survival benefits in HCC within 7 years, while the tumor purity showed a reverse trend. DEGs based on both immune and stromal scores primarily affected the cytokine-cytokine receptor interaction signaling pathway. Among the DEGs, three genes (CASKIN1, EMR3, and GBP5) were found most significantly associated with survival. Moreover, the expression levels of CASKIN1, EMR3, and GBP5 genes were significantly correlated with immune/stromal/Estimate scores or tumor purity and multiple immune cell infiltration. Among them, GBP5 genes were highly related to immune infiltration.
Conclusion: This study identified three key genes which were related to the TME and had prognostic significance in HCC, which may be promising markers for predicting immunotherapy outcomes.

Introduction

Hepatocellular carcinoma (HCC), one of the digestive tract cancers, is also the most common primary liver cancer (1). Hepatocarcinogenesis is a multistep and complex biological process in which many signaling cascades are altered, resulting in heterogeneous molecular profiles and ultimately in tumorigenesis, progression, and metastasis (2). Surgical treatment and chemotherapy are the main therapies for HCC (3), but incidence rates of HCC are continuing to grow, and the probability is rising faster than any other cancer in both men and women (4). A systematic analysis for the Global Burden of Disease Study (GBD) has shown that the incidence and mortality of HCC rank among the top 10 cancers, and death in adults with cirrhosis is the leading cause for the mortality of HCC (5). With the improvement of medical standards, the treatment of HCC has indeed made progress. But currently, HCC treatment is still a global research hotspot, and more and more attention has been paid to cancer immunotherapy, one of the most promising methods for cancer treatment (6). Besides, studies have shown that immune tolerance and escape in the immunosuppressive microenvironment of HCC can be promoted by multiple mechanisms (7, 8). Therefore, it is essential to understand the microenvironment of HCC.

The tumor microenvironment (TME), which is comprised of a mixture of immune cells, stromal cells, cancer cells, the intricate cytokine and chemokines environment, and other components (6, 9), is a dynamic system. Immune cells and stromal cells within the TME are the two main types of non-tumor components which are considered to play important role in the diagnosis and prognosis of tumors (10). Evidence from studies indicates that stromal cells within the TME are genetically stable and are attractive therapeutic targets with reduced risk of resistance and tumor recurrence (11). In addition, due to the dysregulation of the metabolic activity of tumor cells, tumor-infiltrating immune cells usually experience metabolic stress, resulting in an impaired anti-tumor immune response (12). A multi-target approach that simultaneously suppresses TME components may provide a more effective method of treating cancer (13). Therefore, understanding the TME is critical for inhibiting tumorigenesis, invasion, and metastasis, and in effectively managing the immune response (14–16).

Bioinformatics resolve the problems of biology through the methods of applied mathematics, informatics, statistics, and computer science (17, 18). At the same time, with the growth of the amount of biological tumor data, bioinformatics is essential for the storage, analysis, and visualization of cancer immunotherapy data (19, 20). Its rapid development has provided a user-friendly and convenient platform for researchers, guiding the implementation of basic experiments (21, 22). In 2013, Yoshihara et al. created a method to infer the ratio of stromal cells and immune cells in malignant tumors through gene expression signatures that can be derived from The Cancer Genome Atlas (TCGA)–ESTIMATE algorithm. In addition, the algorithm can also predict tumor purity, which helps understand the influence of the microenvironment on neoplastic cells (23). Researchers have extensively verified and confirmed the accuracy of the prediction. In recent years, the ESTIMATE algorithm has been applied to glioblastoma (10), gastric cancer (16), clear cell renal cell carcinoma (9), colon cancer (24), and so on. However, the application of HCC remains to be elucidated.

In this study, we obtained transcriptome data of HCC from the TCGA database, and analyzed the immune/stromal/Estimate scores and tumor purity within the microenvironment using the ESTIMATE algorithm. The relationship between immune/stromal/Estimate scores and tumor purity with survival and clinical parameters were explored. Hub genes associated with immune and stromal scores and with prognostic values were chosen.

Materials and Methods

Data Processing

RNA sequencing (RNA-Seq) and clinicopathological data of patients with HCC were downloaded from Genomic Data Commons

(GDC) database (https://portal.gdc.cancer.gov/) (25, 26). About 374 tumor samples were used for the analysis. ESTIMATE algorithm was used to calculate the stromal/immune/Estimate scores and tumor purity using “estimate” package through the R project (http://r-forge.r-project.org; repos=rforge, dependencies=TRUE) (23). The Human Protein Atlas (HPA) (https://www.proteinatlas.org/about/download) was used to verify the immunohistochemistry staining of genes. Kaplan–Meier plotter (http://kmplot.com/analysis/) was applied to assess the prognostic value of the biomarkers.

Differential Expression Analysis

The median was set as the cut-off value of the immune and stromal scores. Patient samples were divided into two groups, respectively, namely the high-score group and the low-score group. Differential expression analysis was analyzed on the matrix of the sample using the R package, LIMMA. The filtering conditions for the differential genes were as follows: Fold Change > | ±1|, with adjusted p < 0.05.

Enrichment Analysis

The differentially expressed genes (DEGs) from differential expression analysis that meet the conditions were used for enrichment analysis, including Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analysis via Database for Annotation, Visualization, and Integrated Discovery database (DAVID) function annotation tool (https://david.ncifcrf.gov) (27). GO contains biological processes (BPs), cell components (CCs), and molecular functions (MFs), and the signaling pathways were identified by considering both p-value and count number.

Prognostic Gene Selection by Cox Regression Analysis and Least Absolute Shrinkage and Selection Operator (LASSO)

In order to select the survival-related genes, univariate Cox regression analysis was performed on DEGs. At the same time, in order to prevent overfitting and increase the credibility of selecting core genes, we used the least absolute shrinkage and selection operator (LASSO) Cox regression model for signature construction (28, 29). LASSO regression modeling was conducted using the R package, glmnet (30). Finally, multivariable Cox regression analyses were used for feature selection, and calculating hazard ratios (HRs) with 95% confidence intervals (CIs) (31). Genes with p < 0.05 was chosen for further analysis.

Tumor Immune Infiltration Through RNA-Sequencing Expression Profiling Data

Single sample gene set enrichment analysis (ssGSEA) algorithm (https://doi.org/10.1016/j.celrep.2016.12.019, https://doi.org/10.21203/rs.3.rs-33230/v1) is scored based on 29 published immune-related genes (the immune gene set includes immune cell types, functions, and pathways) to quantify the immune infiltration level of 29 immune signatures in each HCC sample using the R package, GSVA (32), and the scores were standardized for each individual immune cell type. Besides, Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) was used to analyze the correlation between six kinds of tumor-infiltrating immune cells and selected hub genes.

Statistical Analysis

Student's t-test is analyzed for comparison between two groups and one-way ANOVA is used to compare multiple groups. Survival analysis was performed using Kaplan–Meier curve by using the R packages, such as survival and survminer, and the p-value was calculated using the log-rank test. Chi-square test was applied to test the association of the expression of the three hub genes with clinicopathological parameters. The value of p < 0.05 was considered statistically significant.

Results

Stromal/Immune/Estimate Scores and Tumor Purity Were Significantly Associated With Prognosis and Tumor Mutation Burden

Transcriptome data and clinical information of HCC were downloaded from the TCGA database and integrated. Totally, 369 cancer samples were used for survival analysis (Supplementary Table 1). At the same time, in order to better understand the impact of tumor infiltrating immune and stromal cells on prognosis, we calculated the immune/stromal/Estimate scores and tumor purity based on the ESTIMATE algorithm, which helps to quantify the immune and stromal components in HCC. Sample information of overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI), and progression-free interval (PFI) were collected (Supplementary Table 2). The results in Figure 1A showed that the four survival periods were significantly correlated with stromal/immune/Estimate scores and tumor purity. Elevated stromal, immune, and Estimate scores were significantly correlated with better prognosis, while higher tumor purity was significantly associated with poor patient survival rate. In addition, we also combined the clinicopathological parameters of HCC to analyze the correlation with immune score, stromal score, Estimate score, and tumor purity (Supplementary Table 3). As can be seen from Supplementary Figure 1, patients with metastasis tended to have higher immune and Estimate scores and lower tumor purity. For tumor and pathological stage, patients with increased stage tended to have high stromal, immune and Estimate score and lower tumor purity.

Figure 1

Analysis of survival and Tumor mutation burden (TMB) associated with immune/stromal/Estimate scores and tumor purity. (A) Kaplan–Meier survival analysis based on immune/stromal/Estimate/scores and tumor purity. Overall survival (OS), disease-specific survival (DSS), disease-free interval (DFI), progression-free interval (PFI). (B) The relationship of immune/stromal/Estimate/scores and tumor purity with tumor mutation burden (TMB). The samples were divided into high and low groups according to the median of score. *p < 0.05, **p < 0.01, and ***p < 0.001 between the two groups.

Tumor mutation burden (TMB) is the total number of somatic gene coding errors, base substitutions, gene insertion, or deletion errors detected per million bases (33). It is a quantitative biomarker that reflects the total number of mutations carried by tumor cells (34). In order to explore the relationship between TMB and immune/stromal/Estimate scores, tumor purity, and the mutation data of HCC from the TCGA database was downloaded for calculating the TMB. Then we divided the tumor samples into high and low groups according to the median of immune/stromal/Estimate scores and tumor purity, respectively. As can be seen from Figure 1B, TMB was significantly higher in the immune/stromal/Estimate scores in the low group, whereas it was higher in tumor purity in the high group.

DEGs Based on Immune and Stromal Scores and Their Associated Pathways Were Identified

To find out novel genes in HCC microenvironment associated with both immune and stromal scores, we performed RNA sequencing (RNA-Seq) differential expression analysis of 374 HCC cases from TCGA cohort. First, we grouped the samples based on the median of the immune score and the stromal score, and then we conducted a difference analysis between the high and low group samples. The DEGs were displayed in the heatmap and volcano plot in Figures 2A,B. The results showed that there were 1,065 significantly upregulated genes and 105 significantly downregulated genes between high and low immune scores. Based on the difference analysis between the high and low stromal scores, 1,597 significantly upregulated genes and 112 significantly downregulated genes were obtained. By overlapping the DEGs of immune and stromal scores, we obtained a total of 896 genes (850 upregulated genes and 46 downregulated genes) for functional enrichment analysis including GO and KEGG (pathways; Figure 2C). The top 10 functional annotations of GO analysis are shown in Figure 2D. DEGs were mostly enriched in receptor activity, plasma membrane, and immune response. In addition, the enriched pathways of DEGs are displayed in Figure 2D. Considering the number of genes enriched in each pathway and the FDR value, the most important pathway was cytokine–cytokine receptor interaction (FDR = 4.24E-19). The interaction of DEGs enriched in the cytokine–cytokine receptor interaction signaling pathway is shown in Figure 3.

Figure 2

Differentially expressed genes (DEGs) based on immune and stromal scores of tumor microenvironment (TME) and their functional annotations in HCC. (A) Heatmaps and volcano plot of the DEGs of stromal scores (p < 0.05, fold change > |±1|). (B) Heatmaps and volcano plot of the DEGs of immune scores (p < 0.5, fold change > |±1|). (C) 46 common downregulated genes and 850 common upregulated genes of both stromal and immune scores were shown by a Venn diagram, and a total of 896 significantly different genes were obtained. (D) The selected DEGs were used for Gene Ontology (GO)-enrichment analysis, biological process (BP), cellular component (CC), and molecular function (MF). Top 10 GO terms were displayed, respectively. (E) The selected DEGs were used for Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis via Kyoto Encyclopedia of Genes and Genomes (DAVID). Considering both the p-value and count number, the optimal pathway was determined.

Figure 3

Pathway diagram showing the interaction of DEGs in the cytokine–cytokine receptor interaction pathway. Alteration frequencies of each gene were represented by the color intensity.

Identification and Validation of Optimal Prognostic Biomarkers in HCC

Through a well-known mathematical model, the prognostic DEGs related to the stromal–immune score that can be used as independent prognostic factors in patients with HCC was identified. We performed a univariate Cox regression analysis on the 896 DEGs related to both stromal and immune scores. Eighty-nine genes with a p < 0.5 were included for further analysis (Supplementary Table 4). To avoid overfitting the variables, 18 genes were screened using LASSO regression analysis (Figure 4A). Finally, three genes including guanylate binding protein 5 (GBP5), adhesion G protein-coupled receptor E3 (EMR3) and CASK interacting protein 1 (CASKIN1) were found significantly associated with HCC prognosis by multivariate Cox regression analysis (Figure 4B).

Figure 4

Screening and verification of prognostic genes in HCC. (A) Least absolute shrinkage and selection operator (LASSO) regression analysis was used to screen for genetic variables. The dotted line indicates the number of genes after screening. (B) Multivariate Cox regression analysis was used to further screen genes that can be used as independent prognostic factors. The value, p < 0.5 was considered statistically significant. (C) Survival verification of the three selected genes by Kaplan–Meier plotter database in HCC. Prognostic indicators include overall survival (OS), progression-free survival (PFS), relapse-free survival (RFS), disease-specific survival (DSS).

We then further verified the prognostic value of the three genes in Kaplan– Meier plotter. Survival analysis of OS, PFS, relapse-free survival (RFS), and DSS according to the expression level of the three genes were performed (Figure 4C). Our results showed that all three genes were significantly associated with the four survival parameters. High expression of the three genes predicted prolonged PFS, RFS, and DSS. High expression of EMR3 was also associated with better OS, while high expression of GBP5 was associated with better OS within 80 months and high expression of CASKIN1 was associated with worse OS before 80 months. To further validate the results, we downloaded and sorted out the GSE76427 gene expression and clinical data of a cohort of 115 HCC cases and 52 adjacent non-tumor tissue from the gene expression omnibus (GEO) database. In addition, we also collected the LIRI–JP dataset in the International Cancer Genome Consortium (ICGC) database, including 243 HCC samples and 202 adjacent samples. Similar to TCGA results, the expression of CASKIN1 gene in HCC is significantly higher than that of normal tissues, and the high expression of GBP5 gene is associated with good OS in the GEO and ICGC data analysis (Supplementary Figure 2).

The Expression Level of GBP5, EMR3, and CASKIN1 and Their Association With Clinicopathological Parameters in HCC

To further confirm the importance of GBP5, EMR3, and CASKIN1 in HCC, Human Protein Atlas (HPA) database was used to compare their protein expression in normal and HCC tissues. As demonstrated in Figure 5A, GBP5 and CASKIN1 were highly expressed in HCC tissue, while EMR3 was downregulated in HCC. At the same time, we also used TCGA database to compare their expression level. CASKIN1 and GBP5 mRNA expression level was significantly increased in HCC tissue compared to adjacent normal tissues, while significantly decreased EMR3 mRNA expression was found in HCC samples compared to normal samples (Figure 5B). Moreover, we also analyzed the association of the three genes with clinicopathological parameters of HCC according to the median of their expressions (Figure 5C). Low expression of EMR3 was often found in low tumor and pathological grade and there was increased proportion of its high expression in late stage. Similar result was observed for pathological stage and tumor grade for CASKIN1 expression.

Figure 5

The expression level of three genes in HCC. (A) The immunohistochemistry (IHC) results from the Human Protein Atlas (HPA) was used to detect the protein level of three genes in normal and tumor tissues. (B) Comparison of the expression levels of CASKIN1, EMR3, and GBP5 genes in HCC tissues and adjacent normal tissues form the Cancer Genome Atlas (TCGA) database. The values of *p < 0.5, **p < 0.01, and ***p < 0.001 between the two groups. (C) Chi-square test of the clinical parameters according to the median of the expressions of the three genes.

The Association of GBP5, EMR3, and CASKIN1 Expression With Immune Cell Infiltration

After identifying the prognostic value and expression level of CASKIN1, EMR3, and GBP5, we performed correlation analysis between CASKIN1, EMR3, and GBP5 expression levels and immune/stromal/Estimate scores, and tumor purity in HCC, respectively. As shown in Figure 6A, CASKIN1 expression was negatively correlated with the immune/stromal/Estimate scores and positively correlated with tumor purity. In contrast, the expression of both EMR3 and GBP5 were positively correlated with immune/stromal/Estimate scores, while negatively correlated with tumor purity. Moreover, the association of GBP5 expression with immune/stromal/Estimate scores and tumor purity was stronger than the other two genes.

Figure 6

Immune infiltration related to CASKIN1, EMR3, and GBP5. (A) Correlation analysis between CASKIN1, EMR3, and GBP5 mRNA expression levels and immune/stromal/Estimate scores and tumor purity in hepatocellular carcinoma (HCC). (B) The relationship between the CASKIN1, EMR3, and GBP5 gene expression and the infiltration level of six types of immune cells in HCC via Tumor Immune Estimation Resource (TIMER) database. Partial Spearman's correlation and statistical analysis were performed. (C) Single sample gene set enrichment analysis (ssGSEA) algorithm was used to obtain the immune-infiltration levels of 29 immune cells. The correlation between the expression levels of CASKIN1, EMR3, and GBP5 expression levels and the infiltration levels of 29 immune cells was displayed using Lollipop Chart.

To study the relationship between the expression of the three genes and immune infiltration level for HCC, scatter plots are shown with partial Spearman's correlation and statistical significance in Figure 6B. The expression of GBP5 was significantly associated with purity (R = −0.461). In addition, elevated EMR3 and GBP5 were significantly correlated with B cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell infiltration (p < 0.05) and a general increase in the immune infiltration level (R > 0.3). It is worth noting that there is a significant positive correlation between the infiltration levels of six kinds of immune cells and the GBP5 gene expression level. Then, based on the gene expression data, we enriched the proportion of 29 immune cells in each patient with HCC through the ssGSEA algorithm, and finally obtained the infiltration level of various immune cells. The correlation of the infiltration level of the immune cells with the expression levels of three genes was calculated. As shown in Figure 6C, the CASKIN1 gene was mainly negatively correlated with immune cell infiltration and the GBP5 gene showed a positive correlation with 29 immune cells and was highly correlated with 12 immune cells (R > 0.3), including T-cell co-inhibition, CD8+_T cells, check-point, inflammation-promoting, TIL, Th1 cells, APC co-inhibition, cytolytic activity, HLA, pDCs, Treg, and Th2 cells. EMR3 gene showed both positive and negative correlation with 29 immune cells.

Discussion

A large number of studies have shown that the TME plays a significant role in the occurrence, development, and metastasis of tumors (11, 13, 35). With the rapid development of bioinformatics based on tumor immunotherapy and microarray sequencing, researchers are increasingly using statistical algorithms to explore new targets for immunotherapy of HCC (19, 36, 37), including ESTIMATE algorithm, and some progress has been made (38–41).

In this study, the ESTIMATE algorithm was used to obtain the immune/stromal/Estimate scores and tumor purity of TME in (HCC). To explore the impact of the immune/stromal/Estimate scores and tumor purity on survival, we collected four types of survival data: OS, DSS, DPI, and PFI. The results showed that the differential immune/stromal/Estimate scores and tumor purity scores significantly affect survival rates. Remarkably, the high score groups of stromal/immune/Estimate scores were significantly associated with longer OS in patients with HCC among the four survival types within 7 years (Figure 1A). These findings are consistent with previous studies showing that the immune/stromal scores were significantly related to OS (40, 42). Findings in multiple cancer types revealed that TMB may play an important role in tumor immunotherapy (43–45), including bladder cancer, colorectal cancer, and non-small cell lung cancer. Therefore, we expect that the immune/stromal/Estimate scores and tumor purity in the microenvironment of HCC are related to TMB. We divided the samples into high and low groups according to the median of the TMB value to compare the correlation between the immune/stromal/Estimate scores, tumor purity, and TMB. The results showed that the higher the immune/stromal scores, the lower is the TMB value. However, the result for tumor purity was the opposite (Figure 1B). This means that the more immune/stromal cells in HCC, the harder it is to identify cancer cells (46). Previous studies have revealed that high TMB predicted worse patient outcomes than those with low TMB in patients with HCC. This is consistent with our survival results (47). In addition, we combined the immune/stromal/Estimate scores and tumor purity with the clinical parameters of HCC, such as metastasis, OS, and grade. However, there was no significant difference (Supplementary Figure 1).

Next, we grouped the samples according to the median of immune score and stromal score to find DEGs. Based on the stromal score, we got 1,579 significantly upregulated and 112 downregulated genes, and based on immune score, 1,046 significantly upregulated genes and 105 significantly downregulated genes were found (Figures 2A,B). We collected 46 DEGs that were jointly downregulated and 850 DEGs that were jointly upregulated in both immune and stromal scores for functional enrichment analysis (Figure 2C), including GO functional annotation analysis (Figure 2D) and KEGG pathway enrichment analysis (Figure 2E). The results indicated that these genes were mainly enriched in receptor activity, plasma membrane, and immune response. It can be seen from the biological process (BP), that the DEGs were closely related to immune response. Consistent with previous studies (48–50), this evidence proved that TME plays a vital role in the immunotherapy of HCC. At the same time, these genes were significantly enriched in the cytokine–cytokine receptor interaction pathway. A simplified pathway diagram in Figure 3 showed DEGs enriched in the pathway.

To evaluate the prognostic significance of these DEGs in HCC, univariate Cox regression analysis, LASSO regression analysis (Figure 4A), multivariate Cox regression analysis (Figure 4B), and Kaplan–Meier survival analysis (Figure 4C) were performed. The results demonstrated that CASKIN1, EMR3, and GBP5 were the most significant prognostic markers. Previous study has demonstrated that EMR3 is one of the adhesion G protein-coupled receptors (aGPCRs), which can be used as a modulator of immune cell function (51). It has prognostic significance in Dukes'B colon cancer (52) and glioblastoma (53). Besides, the GBP5 gene has been studied in a variety of cancers, including gastric adenocarcinoma (54), skin cutaneous melanoma (55, 56), pancreatic adenocarcinoma (57), and HCC, and GBP5 was one of the key genes in the malignant transformation induced by microcystin-LR (MC-LR) in the cell of HCC (58). In addition, GBP5 promotes immunity in mammals. It also plays an important role in regulating human macrophage pyroptosis and uniquely regulates the induction of apoptosis (59). However, the involvement of CASKIN1 gene in cancer has rarely been studied.

Besides, we tested whether these three genes are abnormally expressed in HCC. The results indicated that the expression of CASKIN1, EMR3, and GBP5 showed significant difference between the tumor and normal samples in both TCGA database and HPA (Figures 5A,B). Then, the relationship of the expression of the three genes and their clinical parameters were analyzed via the Chi-square test (Figure 5C). Among them, the expression of EMR3 showed significant difference in the clinical parameters of tumor and pathological stage. There was also significant difference for CASKIN1 in the pathological stage and grade. These results indicated that CASKIN1 and EMR3 were possibly involved in the progression of HCC. At the same time, we also verified the significance of the three prognostic genes in the GEO and ICGC database (Supplementary Figure 2).

Finally, to further analyze the significance of the three genes in tumor immune infiltration, we analyzed the correlation between the three genes and the immune/stromal/Estimate scores and tumor purity, respectively (Figure 6A). At the same time, the TIMER database was applied to assess the correlation between the expression of the three genes and the infiltration scores of six immune cell types. In both the analyses, the GBP5 gene showed a strong correlation with the degree of immune infiltration (Figure 6B). Next, the ssGSEA algorithm was used to evaluate RNA-Seq expression profile data to detect the infiltration of immune cells in tumor tissues of HCC. The correlation analysis between the expression of CASKIN1, EMR3, and GBP5 genes and the immune infiltration scores of 29 immune cell types was performed (Figure 6C). The results showed that GBP5 gene expression exhibited strong positive correlation with 12 kinds of immune cells (R > 0.3, p < 0.5), which verified the results in Figures 6A,B, suggesting that GBP5 may be an important target for targeted immunotherapy of HCC.

In summary, our study identified three TME-related prognostic markers in HCC. CASKIN1 was overexpressed in tumor and its high expression was associated with poor OS, while high expression of EMR3 and GBP5 were associated with better survival. However, the prognostic value of the three genes warrants further validation by more clinical data. Importantly, the GBP5 gene was highly expressed in HCC and strongly correlated with immune cell infiltration. It holds a great potential as a candidate for targeted immunotherapy of HCC.

Data Availability Statement

RNA sequencing and clinicopathological data of HCC patients were downloaded from GDC database (https://portal.gdc.cancer.gov/) - https://gdc-hub.s3.us-east-1.amazonaws.com/latest/TCGA-LIHC.htseq_fpkm.tsv.gz.

Author Contributions

SX, JL, and JS wrote the paper. YZha, XW, XY, and ML edited the manuscript. PK, CC, FD, QW, and YZhe prepared and adjusted the figures. TY and ZX designed the study, provided funding, and reviewed the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.