The role of hypoxia-related genes in TACE-refractory hepatocellular carcinoma: Exploration of prognosis, immunological characteristics and drug resistance based on onco-multi-OMICS approach

Transcatheter arterial chemoembolization (TACE) is an effective treatment for hepatocellular carcinoma (HCC). During TACE, chemotherapeutic agents are locally infused into the tumor and simultaneously cause hypoxia in tumor cells. Importantly, the poor effect of TACE in some HCC patients has been shown to be related to dysregulated expression of hypoxia-related genes (HRGs). Therefore, we identified 33 HRGs associated with TACE (HRGTs) by differential analysis and characterized the mutational landscape of HRGTs. Among 586 HCC patients, two molecular subtypes reflecting survival status were identified by consistent clustering analysis based on 24 prognosis-associated HRGs. Comparing the transcriptomic difference of the above molecular subtypes, three molecular subtypes that could reflect changes in the immune microenvironment were then identified. Ultimately, four HRGTs (CTSO, MMP1, SPP1, TPX2) were identified based on machine learning approachs. Importantly, risk assessment can be performed for each patient by these genes. Based on the parameters of the risk model, we determined that high-risk patients have a more active immune microenvironment, indicating “hot tumor” status. And the Tumor Immune Dysfunction and Exclusion (TIDE), the Cancer Immunome Atlas (TCIA), and Genome of Drug Sensitivity in Cancer (GDSC) databases further demonstrated that high-risk patients have a positive response to immunotherapy and have lower IC50 values for drugs targeting cell cycle, PI3K/mTOR, WNT, and RTK related signaling pathways. Finally, single-cell level analysis revealed significant overexpression of CTSO, MMP1, SPP1, and TPX2 in malignant cell after PD-L1/CTLA-4 treatment. In conclusion, Onco-Multi-OMICS analysis showed that HRGs are potential biomarkers for patients with refractory TACE, and it provides a novel immunological perspective for developing personalized therapies.

Introduction

Hepatocellular carcinoma (HCC) accounts for 85% of all liver cancers (Villanueva., 2019). Despite advances in treatment strategies for HCC, the overall 5-years survival rate for patients with HCC remains below 20% (Zheng et al., 2018). Transcatheter arterial chemoembolization (TACE) is a therapy in which drugs are injected into the arteries supplying HCC tissue (Chang et al., 2020). Some studies suggest that TACE-refractory may lead to poor prognosis in HCC patients. It has been shown that TACE procedures can exacerbate hypoxic states (Liu et al., 2016). However, we still lack a multi-omics data-based perspective on the immunological characteristics of hypoxia-associated gene sets in TACE-refractory patients.

Hypoxia is an intrinsic feature of solid tumors due to the imbalance between tumor cell proliferation rate and vascular nutrient supply (Gray et al., 1953). Previous studies have shown that hypoxia can regulate tumor immune microenvironment (TME), such as promoting the recruitment of innate immune cells, and interfering with the differentiation and function of adaptive immune cells (Feng et al., 2022a). Certain cytokines secreted by malignant tumors, especially in hypoxia condition, may induce angiogenesis and metastasis (Zarogoulidis et al., 2014). A retrospective clinical study has shown that high pre-treatment IL-8 levels are a significant predictor of shorter survival and increased refractoriness of TACE (Kim et al., 2015). Therefore, further studies are needed to investigate the hypoxia-related genes (HRGs) that contribute to TACE refractoriness. Importantly, exploring the relationship between the above genes and drug resistance can lead to the development of new therapeutic strategies.

Nowadays, the study of molecular mechanisms based on Onco-Multi-OMICS approach has become one of the most important tools (Feng et al., 2022b; Zhu et al., 2022). Therefore, we searched for hub HRGs contributing to TACE refractoriness and searched for optimal biological markers by combining transcriptome, single cellome, immunome, and whole-exome. Our study also illustrated the immunological characteristics in different risk group and explored their impact on the response to chemotherapy and immunotherapy. In conclusion, our results were beneficial for the management and treatment of TACE-refractory patients.

Materials and methods

Data collection and pre-processing

The mRNA expression profile data of HCC patients were retrieved from TCGA and GEO databases, and the exclusion criteria was as follows: lack of complete follow-up information, 0 days of survival, and repeated sequencing. Supplementary Table S1 showed treatment details for patients in the GSE14520 cohort before exclusion. Finally, 365 tumor samples from the TCGA-LIHC cohort and 221 tumor samples from the GSE14520 cohort were included. Moreover, to study TACE response, we obtained gene expression profile data from GSE104580, which included 100 TACE-responsive samples and 100 TACE-refractory samples (He et al., 2022). Both somatic mutation data and CNV data were obtained from the TCGA-LIHC cohort, including 371 tumor samples. Notably, ‘maftools’ package was used to present the mutation status of each gene. We removed the batch effect between RNA-seq and microarray data by using the ‘sva’ package and made the newly generated gene matrix based on two cohorts as a meta cohort.

Identification of hypoxia-related genes in TACE refractoriness

Differentially expressed genes (DEGs) between different response states were identified by using the ‘limma’ package in the GSE104580 cohort, p < 0.05, with |logFC| > 1 as the threshold. In addition, 1,694 HRGs were extracted from the previous study (Zhang et al., 2020). The above DEGs and HRGs were overlaid to identify the HRGs associated with TACE (HRGTs).

Enrichment analysis

Differential genes between subtypes were analyzed using the ‘limma’ package (adj. p < 0.05, |logFC| > 1). Biological pathways were annotated by using the ‘clusterProfiler’ package for Gene Ontology (GO), Kyoto Gene and Genome Encyclopedia (KEGG). p-value < 0.05 and q-value < 0.05 were considered as significant enrichment pathways. Differences in biological pathways between subtypes were assessed by using ‘GSVA’ algorithm. And KEGG geneset (c2. cp. kegg. v7. 0. symbols. Gmt) was used as the reference gene set with FDR <0.05 as the threshold.

Consistent clustering analysis

In the meta cohort, the prognostic value of each HRGTs was determined by using univariate cox regression analysis. Consensus clustering, an unsupervised clustering method, is a common method to classify subtypes based on the CDF slope was smallest. Consistent cluster analysis and principal component analysis (PCA) were performed to determine whether each subtype was relatively independent of the other subtypes based on prognostic HRGTs (p < 0.05) and prognostic DEGs (p < 0.05). The number of clusters was determined by using ‘conensusClusterPlus’ package. A 1000 replicates with pltem = 0.8 were performed to verify the stability of the subtypes. We used Kaplan Meier analysis and log-rank test to assess the overall survival (OS) of HCC patients in different subtypes.

Identification and validation of risk scores

Modeling and validation were performed by TCGA-LIHC cohort and GSE14520 cohort, respectively. The least absolute shrinkage and selection operator (LASSO) (Feng et al., 2022c) model was used to remove redundant genes from HRGTs. Subsequently, multivariate Cox regression analysis was performed to integrate the coefficients and then establish risk score formulas by gene expression values. Univariate and multivariate Cox regression analyses were used to assess the prognostic value of risk scores across the entire dataset and the external validation dataset. Time-dependent subject operating characteristic (ROC) curves were used to compare the predictive accuracy of risk scores with traditional clinicopathological parameters.

Drug sensitivity analysis

Half maximal inhibitory concentration (IC50) was calculated using the ‘prophetic’ package. Relevant drugs targeting cell cycle, PI3K/mTOR, WNT, and RTK pathways were obtained from the Genome of Drug Sensitivity in Cancer (GDSC) database. Charoentong et al. developed a quantitative scoring scheme called the Immunophenotype Score (IPS) to identify the determinants affecting tumor immunogenicity. IPS is a representative gene associated with immunogenicity calculated using z-score, and our meta cohort’s IPS was calculated from the Cancer Immunome Atlas (TCIA) database (Wu et al., 2018). Moreover, Peng et al. designed a computational architecture, Tumor Immune dysfunction and ejection (TIDE) score (Jiang et al., 2018), to integrate the two mechanisms of tumor immune escape. Our meta cohort’s TIDE score was calculated from the TIDE database.

Single-cell analysis

The HCC single cell dataset was obtained from GSE125499, and single cell expression profile with annotated cell information were obtained from the Tumor Immune Single Cell Hub (TISCH) database (Sun et al., 2021). Finally, we compared the expression changes of CTSO, MMP1, SPP1, TPX2 in different cell types.

Immunological analysis

We used different algorithms to estimate the abundance of immune cells in different samples, such as ssGSEA, TIMER, CIBERSORT, QUANTISEQ, MCP-counter, XCELL and EPIC. Then, ESTIMATE algorithm was used to calculate the immune score and stromal score to reflect the TME status.

Statistical analysis

All statistical analyses were performed using the R software (v.4.1.1). Kruskal-wallis test was used for comparison between groups, χ2 test was used for association between covariates, and Kaplan-Meier method was utilized to compare survival differences between groups. More detailed statistical methods for transcriptome data processing are covered in the above section (Ye et al., 2022). p < 0.05 was considered statistically significant.

Results

Landscape of HRGTs in HCC

A total of 274 DEGs were first identified from the GSE104580 cohort (Figures 1A,C) and overlapped with existing HRGs genes in the database. Finally, 33 HRGTs were identified (Figure 1B). The above genes may play a key role in TACE refractoriness. HRGTs were mutated in 34 of 371 samples with a frequency of 9.16%, most of which had a low mutation frequency (Figure 1D). In addition, Copy number variation (CNV) were prevalent in HRGTs. ORG1 focused on copy number amplification, while CNV deletion frequency was common in CDC20 (Supplementary Figure S1A). The location of HRGTs on the chromosome (Figure 1E). GSE14520, TCGA-LIHC were included in a meta cohort using the ‘combat’ algorithm. The network of HRGTs specifically described the combined gene interactions and their prognostic significance for patients (Figure 2A). Cox regression analysis identified 24 HRGTs were indicative of prognostic significance for HCC patients (Figure 2B).

FIGURE 1

FIGURE 1. Landscape of HRGTs in HCC. (A) The heat map showed a total of DEGs were identified from the GSE104580 cohort. (B) The venn plot showed DEGs overlapped with existing HRGs in the database. (C) The volcano plot showed dysregulation status of DEGs. (D) Mutation landscape of HRGTs in 371 samples. (E) The location of HRGTs on the human chromosome.

FIGURE 2

FIGURE 2. Molecular subtypes based on HRGTs. (A) The network of HRGTs described the combined interactions and prognostic significance. (B) A forest plot showed 24 HRGTs were indicative of prognostic significance. (C) The classification was optimal when the k value = 2. (D) Heat map of distribution of clinicopathological characteristics and molecular subtypes based on HRGTs. (E) PCA plot revealed that the two molecular subtypes had a relative discrete features. (F) Kaplan-Meier analysis of overall survival time in different molecular subtypes.

Identification of molecular subtypes based on HRGTs

The classification was optimal when the k value = 2 (Figure 2C). Two different subtypes were finally identified, with 326 cases in subtype A and 260 cases in subtype B. The heat map showed the distribution of the clinical features of the different subtypes, with most genes significantly overexpressed in subtype A (Figure 2D). PCA plot revealed that the two molecular subtypes had a relative discrete features (Figure 2E). Prognostic analysis revealed a significant survival disadvantage in the subtype B (Figure 2F).

Immune microenvironment and biological pathways in molecular subtypes

The ESTIMATE algorithm reveals that subtype A has a higher immune score, while the stromal score was significantly downregulated compared to subtype B (Figure 3A). In addition, ssGSEA analysis demonstrated the TME status in different molecular subtypes. We discovered that subtype A is probably exhibiting hot tumor characteristics. This was due to a significant rise in activated CD4⁺ T cells, which may have a more active TME (Figure 3B). In addition, we showed the expression of HLA as well as ICI mRNA in different subtypes. Interestingly, the subtype A had higher mRNA expression in most HLAs, such as HLA-A, HLA-B, HLA-C, and HLA-DDA (Figure 3C). Similarly, subtype A had higher mRNA expression in most ICIs, such as CTLA4 (Figure 3D). We made a hypothesis that subtype A would benefit more from immunotherapy. To explore the biological behavior between these different subtypes, we performed Gene set variation analysis. Subtype A showed significant enrichment with cell cycle pathways compared to subtype B (Figure 3E). In addition, we performed a differential analysis between the two subtypes. It was found that the major enrichment pathways of the 496 DEGs identified (Figure 4A) may be associated with biological processes related to oxidative stress, extracellular genes and drug metabolism (Figures 4B,C). Finally, we identified three different regulatory subtypes based on the above DEGs (Figure 4D). Among them, subtype B had the worst prognosis, while subtype C had the best prognosis (Figure 4E). And the above 33 HRGTs were differentially expressed in different subtypes (Figure 4F).

FIGURE 3

FIGURE 3. Immune microenvironment and biological pathways in molecular subtypes. (A) The Box plot of TME score in different molecular subtypes. (B) The Box plot of immune cells score based on ssGSEA in different molecular subtypes. The Box plot of mRNA expression of HLA (C) and ICIs (D) in different molecular subtypes. (E) GSVA analysis in different molecular subtypes using KEGG genesets. *p < 0.05, **p < 0.01, ***p < 0.001.

FIGURE 4

FIGURE 4. Molecular subtypes based on DEGs. (A) The Volcano plot showed DEGs in different subtypes. Heat map of unsupervised clustering analysis. (B) The circle plot of GO enrichment analysis based on DEGs. (C) The bubble plot of KEGG enrichment analysis based on DEGs. (D) The classification was optimal when the k value = 3. (E) Kaplan-Meier analysis of overall survival time in different molecular subtypes based on DEGs. (F) The Box plot of mRNA expression of HRGTs in different molecular subtypes. *p < 0.05, **p < 0.01, ***p < 0.001.

Identification of risk score in HCC

TCGA-LIHC was used as a training cohort with overall survival (OS) as the outcome. The LASSO model was used to remove redundant genes (Figures 5A,B). The coefficients of each gene were obtained by multifactorial Cox regression analysis. A final signature containing 4 HRGTs was obtained. The formula of each patient was as follows: riskscore = (-0.1277 × expression level of CTSO) + (0.1995 × expression level of MMP1) + (0.1061 × expression level of SPP1) + (0.2385 × expression level of TPX2). Using the median value of risk scores in the TCGA cohort, we identified two risk groups for HCC patients: high-risk group, and low-risk group in all cohorts. Where MMP1, SPP1, and TPX2 were significantly highly expressed in the high-risk group, while CTSO was significantly highly expressed in the low-risk group (Figure 5C). Among them, the risk status plot and the survival distribution plot demonstrated the poorer prognosis of HCC patients with higher risk score (Figures 5D,E). The results of the sankey plot showed a strong association between risk subtypes and molecular subtypes. And most patients in the subtype A and low risk group were in alive status (Figures 5F–H). PCA also showed genomic heterogeneity between the two risk groups (Figure 5I). Survival analysis and ROC curves indicated (Figures 5J,K) that risk score had a good prognostic value in both the TCGA-LIHC cohort and the GSE14520 cohort, and that survival was suboptimal in patients with both high-risk subtypes. Moreover, we performed correlation analysis between hub genes and m6A methylation regulators, and interestingly, it was positively correlated with most of the regulators except IGFBP1, IGFBP2, and IGF2BP1 (Supplementary Figure S2).

FIGURE 5

FIGURE 5. Identification and validation of risk model. (A) Determination of the number of regulators using LASSO analysis. (B) Forest plot of multivariate Cox regression analysis. (C) The heat map showed four HRGTs expression in different risk group (TCGA cohort). Risk status plot (D) and the survival distribution plot (E) demonstrated the poorer prognosis of HCC patients with higher risk score. (F) Sankey diagram based on different subtypes. (G) Differences in risk scores between the two molecular subtypes based on HRGTs. (H) Differences in risk scores between the three molecular subtypes based on DEGs. PCA plot (I), Kaplan-Meier analysis (J), ROC curve of 1,3,5 years (K) of different risk groups in the TCGA and GEO cohort. *p < 0.05, **p < 0.01, ***p < 0.001.

Association of risk subtypes with immune microenvironment

We simultaneously applied different algorithms such as TIMER, CIBERSORT, QUANTISEQ, MCP-counter, XCELL and EPIC to estimate the immune cell infiltration status in each samples. Correlation analysis showed that as the risk score increased, the infiltration score of killer immune cells, such as CD4+ T and CD8+ T cell, also increased (Figure 6A). Similarly, there were differences in the distribution of immune cells in the different risk groups. The high-risk group had a more active TME (Figure 6B). In HLA and ICI analysis, the corresponding mRNA expression was higher in high-risk subtypes (Figures 6C,D). Based on whole-exome sequencing data, patients with both high- and low-risk subtypes did not show significant differences in Top mutated genes, which were TP53, CTNNB1, and TTN (Supplementary Figures S1B,C). Considering the importance of tumor mutational burden (TMB) for immunotherapy, we performed a statistical analysis of the TMB differences between the high- and low-risk groups. It was demonstrated that high-risk group had higher TMB, which suggested that they might have a better response to immunotherapy (Figure 6E). Importantly, when low-TMB and low-riskscore are combined, patients will have the best survival advantage (Figure 6F). In addition, we validated our risk signature in immunotherapy cohort (IMvigor210), and the results were consistent with the above findings, namely, high-risk patients had poor survival outcomes, and high-risk patients were more likely to achieve complete remission (CR), as shown in Supplementary Figure S3.

FIGURE 6

FIGURE 6. Characteristics of immune microenvironment in different risk groups. (A) The heat map showed correlation between risk score and immune function score. (B) The heat map showed differences in immune function of different risk groups. (C) The box plot showed differences in immune checkpoint mRNA expression of different risk groups. (D) The box plot showed differences in HLA mRNA expression of different risk groups. (E) The box plot showed differences in TMB score of different risk groups. (F) Survival analysis by combining TMB score and risk score. *p < 0.05, **p < 0.01, ***p < 0.001.

Risk subtypes could reflect drug resistance in HCC patients

We predicted the drug sensitivity of HCC patients in the meta cohort by utilizing the ‘pRRophetic’ algorithm and a ridge regression model. The results showed that targeting cell cycle (CGP-60474, GW843682x, BI-2536, and CGP-082996) (Figure 7A), PI3K/mTOR signaling (JW-7-52-1, MK-2206, and A-443654) (Figure 7B), WNT signaling (CHIR-99021) (Figure 7C), and RTK signaling (Sunitinib and PHA-665752) (Figure 7D) were significantly more effective in high-risk group than in low-risk group. The TIDE score showed that the effectiveness of immunotherapy may be better in high-risk patients (Figure 7E). In addition, the IPS results also demonstrated that the high-risk group seems to have more immunogenic phenotypes (Figure 7F).

FIGURE 7

FIGURE 7. Risk subtypes could reflect drug resistance. (A) The box plot of targeting cell cycle drug in different risk groups, including CGP-60474, GW843682x, BI-2536, and CGP-082996. (B) The box plot of targeting PI3K/mTOR signaling drug in different risk groups, including JW-7-52-1, MK-2206, and A-443654. (C) The box plot of targeting WNT signaling drug in different risk groups, including CHIR-99021. (D) The box plot of targeting RTK signaling drug in different risk groups, including Sunitinib and PHA-665752. (E) The box plot of TIDE score in different risk groups. (F) The box plot of IPS in different risk groups.

HRGTs in single cell levels

We annotated the GSE125499 single cell expression profile file based on the TISCH database, and t-SNE plot demonstrated the subpopulation of different cells (Figure 8A). In addition, the violin plot demonstrated the expression of CTSO, MMP1, SPP1, and TPX2 in different cell types, with SPP1 being more highly expressed in hepatic progenitor (Figure 8B). Interestingly, after PD-L1/CTLA-4 treatment, CTSO, MMP1, SPP1, and TPX2 were significantly up-regulated in tumor cells (Figure 8C). Finally, we showed the changes in the proportion of different cell types before and after immunotherapy (Figure 8D). The above data suggest to us that four HRGTs involved in risk signature may have a role in reflecting the response to immunotherapy.

FIGURE 8

FIGURE 8. HRGTs in single cell levels. (A) t-SNE plot demonstrated the subpopulation of different cells. (B) The violin plot demonstrated the expression of CTSO, MMP1, SPP1, and TPX2 in different cell types. (C) The violin plot demonstrated the expression of CTSO, MMP1, SPP1, and TPX2 in different treatment groups. (D) The histogram showed the changes in the proportion of different cell types before and after immunotherapys.

Discussion

Primary liver cancer is one of the sixth most common cancers worldwide and is a common tumor of the digestive system with high aggressiveness and poor prognosis (Choi et al., 2020). Since HCC is not sensitive to conventional radiotherapy, surgery becomes the main treatment method. TACE is the treatment of choice for intermediate-stage hepatocellular carcinoma (Morise et al., 2014). In recent years, however, TACE refractoriness has become a thorny issue and has received increasing attention. This is because TACE accompanied by tumor ischemia plays a dual role in the treatment of HCC. Initially, TACE induces tumor necrosis by blocking the vasculature from the hepatic artery to the HCC. However, TACE also stimulates angiogenesis by inducing hypoxia thereby promoting tumor recurrence and metastasis (Kenji et al., 1997). Tumor angiogenesis and invasiveness by TACE have been found to be mediated by hypoxic signaling, which has been effectively inhibited by antiangiogenic therapies (Liu et al., 2016). However, there are few studies related to HRGs associated with TACE. Although a series of studies have identified predictors or models associated with TACE refractoriness, no studies explore the relevance of HGRTs to the immune microenvironment, prognosis and drug resistance.

Onco-Multi-OMICS approach have been commonly used to discover potential biomarkers (Chakraborty et al., 2018). To date, few studies have constructed prognostic models based on combinations of multiple HRGs in TACE-refractory HCC. Importantly, genetic features and clinical characteristics have performed unsatisfactorily in predicting survival outcomes for TACE-refractory patients. Therefore, we aimed to explore a new HRGTs-based risk stratification and propose potential therapeutic targets for HCC patients. Tumor hypoxia promotes the growth of tumor cells and their transformation to a malignant phenotype. The exploration of hypoxia has opened new perspectives for HCC. Hypoxia is a typical hallmark of TME in almost all solid tumors, caused by rapid and uncontrolled tumor proliferation and inadequate blood supply (Graham and Unger, 2018). Under hypoxic conditions, HIFs bind to transcriptional co-activators and hypoxia response elements to increase the expression of a cascade of target genes, thereby regulating various biological processes, including proliferation, metabolism, angiogenesis, migration and invasion. In addition, hypoxia increases the resistance of tumor cells to chemotherapy, radiotherapy and even immunotherapy (Li et al., 2004). It can inactivate effector cytokine production by inhibiting T cell proliferation and function. Therefore, it is important to fully understand the effects of hypoxia on TACE. In our study, a total of 274 DEGs were first identified in the GSE104580 cohort and overlapped with existing HRGs genes in the database. Finally, identifying 33 HRGTs that may have played a key role in TACE refractoriness. Patients were classified into different subtypes according to the expression of prognostic HRGTs and DEGs. The ESTIMATE algorithm showed that subtype A had a higher immune score, and subtype A had higher mRNA expression of most HLAs and ICIs.

ICIs therapy has been shown to be a highly effective agent for the treatment of HCC(Graham and Unger, 2018). However, it is unclear how to identify those who may benefit most from ICIs therapy. Hypoxia promotes tumor progression in different ways, including proliferation, metabolism, angiogenesis and migration, and improves resistance to ICIs(Bao and Wong, 2021). In addition, many factors, especially in TME, can influence the effectiveness of ICIs(Zhang and Zhang, 2020). The TME is a complex and integral component of cancer, containing tumor cells, stromal cells, inflammatory cells, fibroblasts, metabolites and cytokines. To investigate the value of risk subtypes in TME status and immunotherapy, multiple algorithms were used simultaneously in the immune cell analysis to estimate the immune cell infiltration score in different samples. Correlation analysis showed that as the risk score increased, the infiltration fraction of killer immune cells, such as CD4+ T and CD8+ T cell, also increased (Li et al., 2021). And the high-risk group had more active TME. Corresponding mRNA expression was higher in the high-risk subtype in HLA and ICIs analysis. For drug resistance, our study suggested that the high-risk group may have a better response to immunotherapy. We used the pRRophetic algorithm to predict drug sensitivity of HCC patients in different risk groups. The results showed that drugs targeting cell cycle, PI3K/mTOR signaling, WNT signaling, and RTK signaling were more effective in high-risk patients. Importantly, the IPS results demonstrated that the high-risk group seems to have more immunogenic phenotypes.

For the four HRGTs involved in the risk signature, we found that all of them were associated with tumor immunity. Secretory phosphorylated protein 1 (SPP1) is a secreted multifunctional phosphorylated protein that specifically binds and activates matrix metalloproteinases (MMPs) in cancer (Chen et al., 2019a). Its main biological functions are involved in immune response, biomineralization and tissue remodeling. SPP1 has also been implicated in cell growth, proliferation, migration, apoptosis and chemotaxis. Previous studies have demonstrated that SPP1 is overexpressed in a variety of cancers and can be used to predict chemotherapy prognosis, such as ovarian cancer (Zeng et al., 2018), glioblastoma (Kijewska et al., 2017), HCC (Liu et al., 2022) and gastric cancer (Chen et al., 2018a). MMP1 is a member of a family of zinc-dependent endopeptidases involved in wound healing, inflammation, cancer and angiogenic remodeling of the extracellular matrix (ECM) (Chen et al., 2019b). It has been shown to be closely associated with migration and invasion in many cancers. mmp1 promotes cell cycle acceleration in cancer cells by activating the cdc25a/CDK4-cyclin D1 and p21/cdc2-cyclin B1 complexes (Yu et al., 2021). A newly discovered mechanism of MMP1 in tumor promotion is by activating PAR1 to cleave downstream oncogenic signaling pathways (Huang et al., 2018). This is expected to be a promising strategy to address the TACE refractoriness. TPX2 has been identified as an oncogenic factor in a variety of cancers. For example, upregulated expression of TPX2 enhances breast cancer metastasis by mediating MMP2 and MMP9 expression (Tan et al., 2019). In addition, TPX2 can inhibit cell proliferation and enhance apoptosis by blocking the PI3k/AKT/p21 pathway and activating the p53 pathway in breast cancer (Chen et al., 2018b). It has been shown that TPX2 is highly expressed in HCC tissues. TPX2 expression is associated with the infiltration status of immune cells in HCC involving B, CD4⁺T and CD8⁺ T cells, neutrophils, macrophages and DCs(Zhu et al., 2020). In addition, CDK5-mediated stabilization of TPX2 promotes HCC tumorigenesis (Wang et al., 2019). Clearly, these studies suggest that TPX2 is an unfavorable marker for HCC and holds promise as a therapeutic target for TACE refractoriness. CTSO is a cysteine protease that has been shown to have both extracellular and intracellular functions. This class of proteases mediates intracellular protein catabolism and selectively activates extracellular protein degradation, macrophage function and bone resorption (Shi et al., 1995). The role in cancer therapeutic resistance is an emerging area of interest.

In our study, different hypoxic patterns present different biological processes, signaling pathways and immune features. Based on the parameters of the risk model, we determined that high-risk patients have a more active immune microenvironment, and HRGs are potential biomarkers for TACE-refractory patients. Especially, it may be an independent prognostic factor for HCC patients. However, our study has some limitations. Firstly, we should use advanced artificial intelligence models rather than traditional machine learning models such as Random Forest (RF) or LASSO models. However, for clinical applications, machine learning models with coefficients may be more helpful to clinicians. The clinician can calculate the survival risk of each patient from the mRNA expression and coefficient, however, more advanced deep learning models are a ‘black box’. Moreover, due to the limitation of laboratory conditions, we have no more time to conduct in vivo or vitro experiments, and we will validate the mechanism of four hub genes in TACE-refractory patients in the future. In conclusion, our study will provide a novel immunological perspective for the development of treatment options for TACE-refractory HCC.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

SF and XC conceived and designed the study. SF was responsible for materials. XC, JL, and LF drafted the article. XW and YL revised the article critically. All authors had final approval of the submitted versions.

Funding

This work was partially supported by funding from the National Natural Science Foundation of China (No. 81672944), Natural Science Foundation of Jiangsu Province (Nos. BK20211296 and BK20200843), the Priority Academic Program Development of Jiangsu Higher Education Institutions (Integration of Chinese and Western Medicine), and Innovation Team of Six Talent Peaks Project in Jiangsu Q12 Province (No. WSN-062).

Acknowledgments

We are grateful to Feng from southeast university for bioinformatics analysis, and computer resource from state key laboratory of bioelectronics, school of biological science and medical engineering, southeast university.

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.1011033/full#supplementary-material

Supplementary Table S1 | Treatment details for HCC patients in GSE14520.

References

Bao, M. H., and Wong, C. C. (2021). Hypoxia, metabolic reprogramming, and drug resistance in liver cancer. Cells 10 (7), 1715. PMID: 34359884; PMCID: PMC8304710. doi:10.3390/cells10071715

PubMed Abstract | CrossRef Full Text | Google Scholar

Chakraborty, S., Hosen, M. I., Ahmed, M., and Shekhar, H. U. (2018). Onco-multi-OMICS approach: A new frontier in cancer research. Biomed. Res. Int. 2018, 9836256. PMID: 30402498; PMCID: PMC6192166. doi:10.1155/2018/9836256

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, Y., Jeong, S. W., Young Jang, J., and Jae Kim, Y. (2020). Recent updates of transarterial chemoembolilzation in hepatocellular carcinoma. Int. J. Mol. Sci. 21 (21), 8165. PMID: 33142892; PMCID: PMC7662786. doi:10.3390/ijms21218165

CrossRef Full Text | Google Scholar

Chen, L. Z., He, C. Y., Su, X., Peng, J. L., Chen, D. L., Ye, Z., et al. (2018). SPP1 rs4754 and its epistatic interactions with SPARC polymorphisms in gastric cancer susceptibility. Gene 640, 43–50. Epub 2017 Sep 28. PMID: 28962925. doi:10.1016/j.gene.2017.09.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Zhang, H., Zhang, G., Zhong, A., Ma, Q., Kai, J., et al. (2018). Targeting TPX2 suppresses proliferation and promotes apoptosis via repression of the PI3k/AKT/P21 signaling pathway and activation of p53 pathway in breast cancer. Biochem. Biophys. Res. Commun. 507 (1-4), 74–82. Epub 2018 Nov 16. PMID: 30454896. doi:10.1016/j.bbrc.2018.10.164

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Xiong, D., Ye, L., Yang, H., Mei, S., Wu, J., et al. (2019). SPP1 inhibition improves the cisplatin chemo-sensitivity of cervical cancer cell lines. Cancer Chemother. Pharmacol. 83 (4), 603–613. Epub 2019 Jan 9. PMID: 30627777. doi:10.1007/s00280-018-3759-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Peng, S., Cen, H., Lin, Y., Huang, C., Chen, Y., et al. (2019). MicroRNA hsa-miR-623 directly suppresses MMP1 and attenuates IL-8-induced metastasis in pancreatic cancer. Int. J. Oncol. 55 (1), 142–156. PMID: 31115512; PMCID: PMC6561617. doi:10.3892/ijo.2019.4803

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, J., Lee, D., Shim, J. H., Kim, K. M., Lim, Y. S., Lee, Y. S., et al. (2020). Evaluation of transarterial chemoembolization refractoriness in patients with hepatocellular carcinoma. PLoS One 15 (3), e0229696. PMID: 32130270; PMCID: PMC7055892. doi:10.1371/journal.pone.0229696

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, S., Xia, T., Ge, Y., Zhang, K., Ji, X., Luo, S., et al. (2022). Computed tomography imaging-based radiogenomics analysis reveals hypoxia patterns and immunological characteristics in ovarian cancer. Front. Immunol. 13, 868067. PMID: 35418998; PMCID: PMC8995567. doi:10.3389/fimmu.2022.868067

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, S., Xu, Y., Dai, Z., Yin, H., Zhang, K., and Shen, Y. (2022). Integrative analysis from multicenter studies identifies a WGCNA-derived cancer-associated fibroblast signature for ovarian cancer. Front. Immunol. 13, 951582. PMID: 35874760; PMCID: PMC9304893. doi:10.3389/fimmu.2022.951582

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, S., Yin, H., Zhang, K., Shan, M., Ji, X., Luo, S., et al. (2022). Integrated clinical characteristics and omics analysis identifies a ferroptosis and iron-metabolism-related lncRNA signature for predicting prognosis and therapeutic responses in ovarian cancer. J. Ovarian Res. 15 (1), 10. PMID: 35057848; PMCID: PMC8772079. doi:10.1186/s13048-022-00944-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Graham, K., and Unger, E. (2018). Overcoming tumor hypoxia as a barrier to radiotherapy, chemotherapy and immunotherapy in cancer treatment. Int. J. Nanomedicine 13, 6049–6058. PMID: 30323592; PMCID: PMC6177375. doi:10.2147/IJN.S140462

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, L., Conger, A., Ebert, M., Hornsey, S., and Scott, O. (1953). The concentration of oxygen dissolved in tissues at the time of irradiation as a factor in radiotherapy. Br. J. Radiol. 26 (312), 638–648. doi:10.1259/0007-1285-26-312-638

PubMed Abstract | CrossRef Full Text | Google Scholar

He, Q., Yang, J., and Jin, Y. (2022). Development and validation of TACE refractoriness-related diagnostic and prognostic scores and characterization of tumor microenvironment infiltration in hepatocellular carcinoma. Front. Immunol. 13, 869993. PMID: 35493518; PMCID: PMC9043752. doi:10.3389/fimmu.2022.869993

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Li, Y., Guo, Y., Zhang, Z., Lian, G., Chen, Y., et al. (2018). MMP1/PAR1/SP/NK1R paracrine loop modulates early perineural invasion of pancreatic cancer cells. Theranostics 8 (11), 3074–3086. PMID: 29896303; PMCID: PMC5996366. doi:10.7150/thno.24281

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. Epub 2018 Aug 20. PMID: 30127393; PMCID: PMC6487502. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kenji, J., Hyodo, I., Tanimizu, M., Tanada, M., Nishikawa, Y., Hosokawa, Y., et al. (1997). Total necrosis of hepatocellular carcinoma with a combination therapy of arterial infusion of chemotherapeutic lipiodol and transcatheter arterial embolization: Report of 14 cases. Semin. Oncol 24 (2), S6-S–71. PMID: 9151920.

PubMed Abstract | Google Scholar

Kijewska, M., Kocyk, M., Kloss, M., Stepniak, K., Korwek, Z., Polakowska, R., et al. (2017). The embryonic type of SPP1 transcriptional regulation is re-activated in glioblastoma. Oncotarget 8 (10), 16340–16355. PMID: 28030801; PMCID: PMC5369967. doi:10.18632/oncotarget.14092

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S. S., Cho, H. J., Won, J. H., Bae, J. I., Kang, D. R., Lee, J. D., et al. (2015). Interleukin-8 level as a prognostic marker in patients with Hepatitis B virus-associated hepatocellular carcinoma treated with transarterial chemoembolization. Cytokine 76 (2), 449–457. Epub 2015 Jul 8. PMID: 26163999. doi:10.1016/j.cyto.2015.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Teixeira, A. F., Zhu, H. J., and Ten Dijke, P. (2021). Cancer associated-fibroblast-derived exosomes in cancer progression. Mol. Cancer 20 (1), 154. PMID: 34852849; PMCID: PMC863844. doi:10.1186/s12943-021-01463-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Feng, G. S., Zheng, C. S., Zhuo, C. K., and Liu, X. (2004). Expression of plasma vascular endothelial growth factor in patients with hepatocellular carcinoma and effect of transcatheter arterial chemoembolization therapy on plasma vascular endothelial growth factor level. World J. Gastroenterol. 10 (19), 2878–2882. PMID: 15334691; PMCID: PMC4572123. doi:10.3748/wjg.v10.i19.2878

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, K., Min, X. L., Peng, J., Yang, K., Yang, L., and Zhang, X. M. (2016). The changes of HIF-1α and VEGF expression after TACE in patients with hepatocellular carcinoma. J. Clin. Med. Res. 8 (4), 297–302. PMID: 26985249; PMCID: PMC4780492. doi:10.14740/jocmr2496w

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Zhang, R., Deng, J., Dai, X., Zhu, X., Fu, Q., et al. (2022). Construction of TME and Identification of crosstalk between malignant cells and macrophages by SPP1 in hepatocellular carcinoma. Cancer Immunol. Immunother. 71 (1), 121–136. Epub 2021 May 24. PMID: 34028567. doi:10.1007/s00262-021-02967-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Morise, Z., Kawabe, N., Tomishige, H., Nagata, H., Kawase, J., Arakawa, S., et al. (2014). Recent advances in the surgical treatment of hepatocellular carcinoma. World J. Gastroenterol. 20 (39), 14381–14392. PMID: 25339825; PMCID: PMC4202367. doi:10.3748/wjg.v20.i39.14381

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, G. P., Chapman, H. A., Bhairi, S. M., DeLeeuw, C., Reddy, V. Y., and Weiss, S. J. (1995). Molecular cloning of human cathepsin O, a novel endoproteinase and homologue of rabbit OC2. FEBS Lett. 357 (2), 129–134. PMID: 7805878. doi:10.1016/0014-5793(94)01349-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, D., Wang, J., Han, Y., Dong, X., Ge, J., Zheng, R., et al. (2021). TISCH: A comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. 49 (D1), D1420–D1430. PMID: 33179754; PMCID: PMC7778907. doi:10.1093/nar/gkaa1020

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, G. Z., Li, M., Tan, X., Shi, M. L., and Mou, K. (2019). MiR-491 suppresses migration and invasion via directly targeting TPX2 in breast cancer. Eur. Rev. Med. Pharmacol. Sci. 23 (22), 9996–10004. doi:10.26355/eurrev_201911_19566

PubMed Abstract | CrossRef Full Text | Google Scholar

Villanueva, A. (2019). Hepatocellular carcinoma. N. Engl. J. Med. 380 (15), 1450–1462. doi:10.1056/NEJMra1713263

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., Zhao, W., Gao, Y., Zhou, J., Li, H., Zhang, G., et al. (2019). CDK5-mediated phosphorylation and stabilization of TPX2 promotes hepatocellular tumorigenesis. J. Exp. Clin. Cancer Res. 38 (1), 286. PMID: 31272499; PMCID: PMC6610961. doi:10.1186/s13046-019-1297-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J., Zhao, W., Zhou, B., Su, Z., Gu, X., Zhou, Z., et al. (2018). TSNAdb: A database for tumor-specific neoantigens from immunogenomics data analysis. Genomics Proteomics Bioinforma. 16 (4), 276–282. Epub 2018 Sep 15. PMID: 30223042; PMCID: PMC6203688. doi:10.1016/j.gpb.2018.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, F., Liang, Y., Cheng, Z., Liu, Y., Hu, J., Li, W., et al. (2022). Immunological characteristics of alternative splicing profiles related to prognosis in bladder cancer. Front. Immunol. 13, 911902. PMID: 35769470; PMCID: PMC9234272. doi:10.3389/fimmu.2022.911902

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Xu, Z., Guo, J., Yang, K., Zheng, J., and Sun, X. (2021). Tumor-associated macrophages (TAMs) depend on MMP1 for their cancer-promoting role. Cell Death Discov. 7 (1), 343. PMID: 34753916; PMCID: PMC8578434. doi:10.1038/s41420-021-00730-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zarogoulidis, P., Katsikogianni, F., Tsiouda, T., Sakkas, A., Katsikogiannis, N., and Zarogoulidis, K. (2014). Interleukin-8 and interleukin-17 for cancer. Cancer Invest. 32 (5), 197–205. Epub 2014 Mar 26. PMID: 24669909. doi:10.3109/07357907.2014.898156

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, B., Zhou, M., Wu, H., and Xiong, Z. (2018). SPP1 promotes ovarian cancer progression via Integrin β1/FAK/AKT signaling pathway. Onco. Targets. Ther. 11, 1333–1343. PMID: 29559792; PMCID: PMC5856063. doi:10.2147/OTT.S154215

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Tang, B., Gao, J., Li, J., Kong, L., and Qin, L. (2020). A hypoxia-related signature for clinically predicting diagnosis, prognosis and immune microenvironment of hepatocellular carcinoma patients. J. Transl. Med. 18 (1), 342. PMID: 32887635; PMCID: PMC7487492. doi:10.1186/s12967-020-02492-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., and Zhang, Z. (2020). The history and advances in cancer immunotherapy: Understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol. Immunol. 17 (8), 807–821. Epub 2020 Jul 1. PMID: 32612154; PMCID: PMC7395159. doi:10.1038/s41423-020-0488-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, R., Qu, C., Zhang, S., Zeng, H., Sun, K., Gu, X., et al. (2018). Liver cancer incidence and mortality in China: Temporal trends and projections to 2030. Chin. J. Cancer Res. 30 (6), 571–579. PMID: 30700925; PMCID: PMC6328503. doi:10.21147/j.issn.1000-9604.2018.06.01

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, H., Liu, J., Feng, J., Zhang, Q., Bian, T., Li, X., et al. (2020). Overexpression of TPX2 predicts poor clinical outcome and is associated with immune infiltration in hepatic cell cancer. Med. Baltim. 99 (49), e23554. PMID: 33285774; PMCID: PMC7717782. doi:10.1097/MD.0000000000023554

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, Y., Feng, S., Song, Z., Wang, Z., and Chen, G. (2022). Identification of immunological characteristics and immune subtypes based on single-sample gene set enrichment analysis algorithm in lower-grade glioma. Front. Genet. 13, 894865. PMID: 35646050; PMCID: PMC9136245. doi:10.3389/fgene.2022.894865

PubMed Abstract | CrossRef Full Text | Google Scholar