Skip to main content
Download PDF

Exploration and validation of a combined Hypoxia and m6A/m5C/m1A regulated gene signature for prognosis prediction of liver cancer

BMC Genomics volume 24, Article number: 776 (2023) Cite this article

Abstract

Background

It is widely acknowledged that hypoxia and m6A/m5C/m1A RNA modifications promote the occurrence and development of tumors by regulating the tumor microenvironment. This study aimed to establish a novel liver cancer risk signature based on hypoxia and m6A/m5C/m1A modifications.

Methods

We collected data from The Cancer Genome Atlas (TCGA-LIHC), the National Omics Data Encyclopedia (NODE-HCC), the International Cancer Genome Consortium (ICGC), and the Gene Expression Omnibus (GEO) databases for our study (GSE59729, GSE41666). Using Cox regression and least absolute shrinkage and selection operator (LASSO) method, we developed a risk signature for liver cancer based on differentially expressed genes related to hypoxia and genes regulated by m6A/m5C/m1A modifications. We stratified patients into high- and low-risk groups and assessed differences between these groups in terms of gene mutations, copy number variations, pathway enrichment, stemness scores, immune infiltration, and predictive capabilities of the model for immunotherapy and chemotherapy efficacy.

Results

Our analysis revealed a significantly correlated between hypoxia and methylation as well as m6A/m5C/m1A RNA methylation. The three-gene prognosis signature (CEP55, DPH2, SMS) combining hypoxia and m6A/m5C/m1A regulated genes exhibited strong predictive performance in TCGA-LIHC, NODE-HCC, and ICGC-LIHC-JP cohorts. The low-risk group demonstrated a significantly better overall survival compared to the high-risk group (p < 0.0001 in TCGA, p = 0.0043 in NODE, p = 0.0015 in ICGC). The area under the curve (AUC) values for survival at 1, 2, and 3 years are all greater than 0.65 in the three cohorts. Univariate and Multivariate Cox regression analyses of the three datasets indicated that the signature could serve as an independent prognostic predictor (p < 0.001 in the three cohorts). The high-risk group exhibited more genome changes and higher homologous recombination deficiency scores and stemness scores. Analysis of immune infiltration and immune activation confirmed that the signature was associated with various immune microenvironment characteristics. Finally, patients in the high-risk group experienced a more favorable response to immunotherapy, and various common chemotherapy drugs.

Conclusion

Our prognostic signature which integrates hypoxia and m6A/m5C/m1A-regulated genes, provides valuable insights for clinical prediction and treatment guidance for liver cancer patients.

Peer Review reports

Introduction

Liver cancer, the second most deadly malignancy worldwide, is characterized by a poor prognosis [1]. Existing diagnostic tools for liver cancer relying on histological and radiological assessments, often lack precision and practicality [2]. Moreover, clinicopathological features have failed to account for the heterogeneity, among patients, resulting in varying treatment outcomes even among individuals with the same TNM stage [3]. Current staging methods also fall short in guide immunotherapy, which has shown promise in treating liver cancer [4, 5]. Hence, there is pressing need for the development of a novel liver cancer signature capable of predicting prognosis and guiding treatment decisions.

The presence of hypoxia within tumor microenvironment (TME) has been attributed to the considerable distance between the vascular system and tumor cells, coupled with the rapid proliferation of tumor cells [6, 7]. As oxygen is becomes increasingly scarce, tumor tissues establish a new dynamic equilibrium, leading to the formation a hypoxia, hypoglycemic and acidic TME that favors tumor growth [8, 9]. Extensive research has established a strong correlation between hypoxia and rapid tumor progression, metastasis, and drug resistance [10,11,12,13]. Notably, liver cancer is one of the most hypoxic malignant tumors, with a median oxygen content as low as 0.8% [14]. The hypoxic TME induces metabolic reprogramming, mediated by mRNA methylation, which subsequently triggers phenotypic changes in immune cells, creating an immunosuppressive TME [15].

Eukaryotic mRNA undergoes methylation modifications including N6-methyladenosine (m6A), 5-methylcytosine (m5C), and N1-methyladenosine (m1A). Research has revealed that genes involved in the regulation of m6A/m5C/m1A modifications play a pivotal role in shaping the TME, thereby promoting tumor progression [16,17,18,19]. Hypoxia-inducible factors (HIFs) orchestrate cellular adaptation to hypoxic environments by participating in multiple regulatory pathways [20, 21]. HIFs activation governs glucose and lactate metabolism, enabling tumors to thrive in hypoxia conditions while also contributing to the formation of an immunosuppressive TME [22,23,24,25]. Studies have demonstrated that m6A-modifying enzymes influence the methylation and expression of HIFs, thus affecting tumor proliferation. For example, in glioblastoma, hypoxia-induced ALKBH5 stabilizes SFPQ on the CXCL8 gene by clearing m6A methylated NEAT1, promoting the expression of CXCL8/IL8 and facilitating immune evasion [26]. Upregulation of HBXIP can enhance the expression of m6A methylase METTL3, maintaining elevated levels of HIF-α and promoting the malignant proliferation of liver cancer [27]. Hypoxia within the TME alters the methylation levels of m6A/m5C/m1A, subsequently affecting downstream adaptive responses such as immune cell function and tumor behavior [28,29,30]. It has been reported that there is a positive feedback loop between hypoxia and m6A/m5C/m1A methylation in driving malignant tumor proliferation [15].

Therefore, in the present study, we analyzed the interplay between hypoxia and m6A/m5C/m1A regulated genes and established a novel risk prognosis signature for liver cancer that integrated both hypoxia and m6A/m5C/m1A modification. To our knowledge, this represents the first prognostic signature that combines hypoxia and mRNA methylation regulatory genes, offering an enriching approach to clinical management and providing valuable guidance for the neoadjuvant treatment of liver cancer.

Materials and methods

Dataset collection and preprocessing

The flow chart outlining our research is illustrated in Fig. S1. We collected RNA-seq, methylation-seq, and clinicopathological data for liver cancer samples from various sources. Specifically, TCGA-LIHC and ICGC-JP were acquired from Xena (https://xenabrowser.net/datapages/), NODE-OEP000321 were downloaded from NODE (https://www.biosino.org/node/), and GSE59729, GSE41666, and GSE144269 were obtained from GEO (https://www.ncbi.nlm.nih.gov/gds). The clinical information of these data sets is listed in Table 1. To facilitate analysis, we converted all transcriptome sequencing data from count to TPM format and retrieved chip data using R package GEOquery (version 2.64.2) [31]. Subsequently, we standardized all gene expression profiles using the normalizeBetweenArrays function from the R package limma (only for chip data) and log2(X + 1) transformation [32]. The TCGA-LIHC dataset was employed as the training set for model development, while NODE-OEP000321 and ICGC-JP served as the validation sets. All analyses were performed using R (version 4.2.2). GSE144269, a transcriptome sequencing dataset for liver cancer, was used to verify the expression levels of risk genes.

Table 1 Summary of clinical information from the datasets used in the analysis
Full size table

Hypoxia genes were extracted from the hallmark gene set in the Molecular Signatures Database (MsigDb, www.gsea-msigdb.org, version 7.0). The regulatory genes (Writer, Reader, and Eraser) of m6A/m5C/m1A were obtained from the literature [33, 34].

Unsupervised clustering of hypoxia and its correlation with methylation

To assess the hypoxia status of liver cancer patients, we applied an unsupervised clustering algorithm by the R package ConsensusClusterPlus (version 1.60.0) (default parameter) based on the hypoxia hallmark gene set [35]. Additionally, we computed the activation scores for hypoxia and m6A/m5C/m1A regulation using the single sample Gene Set Enrichment Analysis (ssGSEA) method from the R package GSVA (version 1.44.2). We analyzed differences in the expression of m6A/m5C/m1A regulator genes and total methylation levels across different hypoxia clusters. Furthermore, we examined the correlation between hypoxia and the activation of m6A/m5C/m1A regulation [36].

Established a Hypoxia-m6A/m5C/m1A-related-score scoring system to evaluate liver cancer cases

To identify differentially expressed genes (DEGs) associated with hypoxia status more accurately, we intersected DEGs obtained from hypoxia clustering (TCGA-LIHC) and differential oxygen treatment (GSE41666 and GSE59729). We then calculated DGEs using the R package limma (|logFoldChange|> 1, p < 0.01) [32]. Subsequently, we analyzed DEGs between tumor and normal samples using the same package (|logFoldChange|> 1, p < 0.01) and identified gene sets significantly related to m6A/m5C/m1A regulator genes (r > 0.5, p < 0.001). The intersection of these two gene sets yielded a candidate gene set for model establishment using joint Cox regression analysis. Initially, we employed the Univariate Cox regression to identify genes related to overall survival (OS) (selecting genes with p values less than 0.05 after adjusting for p values using the Benjamini & Hochberg correction) via the R package Survival (version 3.5–5). Subsequently, we employed the Least Absolute Shrinkage and Selection Operator (LASSO) penalized Cox regression model to mitigate overfitting among candidate genes using the R package glmnet (version 4.1–7) and Survival. Finally, a risk model was established using Multivariate Cox regression, and risk scores for each patient were computed as follows:

$$\mathrm{risk~score}~({\text{patient}})=\sum_{k=1}^{n} \left({\text{coef}} \, \times {\text{exp}}\right)$$

where n means all samples, exp indicates the expression level for each risk gene, and coef indicates their regression coefficients.

Evaluation and verification of the prognostic signature

We assessed the predictive and generalization capabilities of the Hypoxia-m6A/m5C/m1A-related-Score (HMRs) gene signature in both the training and validation sets. Patients were stratified into high- and low-risk groups based on the median risk score. We analyzed survival differences between these groups using Kaplan–Meier (K-M) curves with the R package survminer (version 0.4.9) and evaluated the predictive capability of the risk score through Received Operating characteristic curve (ROC) analysis. Univariate and Multivariate Cox regression analyses were conducted to assess the clinical prognostic independence of the prognostic risk model. Principal Component Analysis (PCA) and t-Distributed Stochastic Neighbor Embedding (t-SNE) were employed to evaluate the classification ability of the HMRs signature.

Relationships between HMRs gene signature and genomic alterations

It is well-established that cancer patients with different prognoses can distinguish variations in gene expression and mutation patterns. Consequently, we investigated the association between risk groups determined by HMRs and genomic alterations. Initially, we performed an analysis of DEGs between high- and low-risk groups using the R package limma (|logFoldChange|> 0.5, p < 0.05). To gain insights into the functional relevance of these DEGs, we conducted pathway enrichment analysis using Metascape [37]. Subsequently, we utilized the R package progeny (version 1.18.0) to assess the activation levels of 14 typical cancer related pathways between the high- and low-risk groups. We then explored the relationship between genomic instability and risk scores in liver cancer from various angles. This included an analysis of the mutation spectrum and somatic total mutation burden (TMB) using the R package maftools (version 2.12.0) and an examination of copy number variation (CNV) using GISTIC2 [38, 39]. Additionally, we investigated differences in fragmented genomic alterations and Homologous Recombination Deficiency (HRD) between the two risk groups [40, 41]. Cancer stem cells, which exhibit characteristics akin to normal stem cells, play a pivotal role in drug resistance, proliferation, and metastasis of tumors [42]. To better understand their involvement, we compared tumor stem cell indices between the high- and low-risk groups [43].

Relationships between HMRs gene signature and immunocyte infiltration

It is widely acknowledged that the immune microenvironment of tumor tissue has a profound impact on tumor progression and prognosis [44]. Therefore, we examined the relationship between the immune microenvironment of liver cancer and the HMRs gene signature from various angles. Initially, we utilized the R package estimate (version 1.0.13) to compute the immune score of the tumor [45]. Subsequently, we employed the ssGSEA method from the R package GSVA (version 1.44.2) to estimate the proportions of 28 typical immune cell types [36]. To comprehensively assess the immune microenvironment of liver cancer, we employed a variety of immunoinfiltration analysis methods, including TIMER, CIBERSORT, XCELL, QUANTISEQ, MCP-counter, EPIC, and CIBERSORT-ABS. All proportions of immune cells were retrieved from TIMER2.0 (http://timer.comp-genomics.org/) [46,47,48,49,50,51,52]. To evaluate associations between immune cells and risk scores, we employed Wilcoxon signed-rank sum tests and Pearson correlation analyses. Finally, we examined the activation of the immune cycle in liver cancer based on the Tracking Tumor Immunophenotype (TIP) tool (http://biocc.hrbmu.edu.cn/TIP/index.jsp) [53].

Estimate of immunotherapeutic and chemotherapy drug response between HMRs groups

To broaden the model’s applicability, we analyzed its potential impact on clinical practice. Initially, we assessed the predictive power of risk scores for immunotherapy outcomes. We conducted this analysis using the PD-1 treatment cohort from IMvigor210 to explore the relationship between HMRs and immunotherapy response [54]. Additionally, we then evaluated the predictive effect of HMRs on immunotherapy using Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/faq/), tumor inflammation score (TIS), and immunophenoscore (IPS, obtained from The Cancer Immunome Atlas Database, https://tcia.at/about) [55,56,57,58]. Furthermore, we compared differences in common immune checkpoint inhibitors (ICI) and the Human Leukocyte Antigen (HLA) family genes among the HMRs groups. Transcatheter arterial chemoembolization (TACE) serve as a primary treatment for unresectable liver cancer. However, its therapeutic efficacy varies due to the heterogeneity of liver cancer. We analyzed the ability of the risk score to predict TACE response based on GSE104580. Finally, we utilized the Genomic of Drug Sensitivity in Cancer (GDSC) and the R package pRRophetic (version 0.5) to predict the half-maximal inhibitory concentration (IC50) of 8 common liver cancer chemotherapeutic drugs. We compared differences in IC50 values between the HMRs groups by using Wilcoxon signed-rank sum tests [59, 60].

Multi-database validation of risk gene expression and function analysis of candidate gene

To validate the findings, we leveraged multiple independent datasets to analyze the differential expression of risk genes between tumor and normal tissues. Initially, we identified the differentially expressed risk genes in the TCGA-LIHC, ICGC-JP, and GSE144269 datasets. Furthermore, we employed the online tool, The University of Alabama at Birmingham CANcer Data Analysis Portal (UALCAN, https://ualcan.path.uab.edu/index.html), to assess differences in risk gene expression at the protein level [61]. Immunohistochemical results for risk genes were obtained from The Human Protein Atlas (THPA) database (https://www.proteinatlas.org/) [62]. Finally, we examined the expression differences of risk genes in paired and unpaired tumor datasets using the online tool TNMplot (https://tnmplot.com/analysis/) [63].

We also conducted functional analyses of risk genes, including pan-cancer survival and differential expression analysis, using Gene Expression Profiling Interactive Analysis (GEPIA2, http://gepia2.cancer-pku.cn/#index) [64]. Additionally, differences in the methylation levels of genes under different clinical classifications were analyzed using UALCAN. To visualize the methylation sites of risk genes, we employed the R package trackViewer (version 1.36.2) to visualize [65]. Finally, we performed single-cell expression analysis of risk genes, initially assessing the expression levels of risk genes across different cell types in multiple liver cancer single-cell datasets using the Tumor Immune Single-cell Hub 2 (TISCH2, http://tisch.comp-genomics.org/home/) [66]. Subsequently, we determined in which cell types the risk genes were significantly expressed using the PanglaoDB database (https://panglaodb.se/index.html) [67]. Finally, the relationship between the expression levels of risk genes at both protein and transcript levels and the cell cycle were explored using the THPA database [62].

Results

Hypoxia-based clustering of liver cancer was significantly correlated with m6A/m5C/m1A methylation regulation

The results of unsupervised clustering showed that hypoxia gene sets could classify liver cancer into two clusters (Fig. 1A, B), with significant differences in survival between these two clusters (Fig. 1C). Pathway activation analysis based on ssGSEA showed that the activation degree of hypoxia in cluster 2 was significantly higher than that in cluster 1 (Fig. 1D). Almost all m6A/m5C/m1A regulator genes were significantly highly expressed in cluster 2 (Fig. 1E), and the average methylation degree in cluster 2 was significantly higher than that in cluster 1 (Fig. 1F). Pathway activation analysis based on ssGSEA showed that the activation degree of m6A/m5C/m1A in cluster 2 was significantly higher than in cluster 1 (Fig. 1G), and there was a significant correlation between hypoxia pathway activation of m6A/m5C/m1A regulation pathway (Fig. 1H, r = 0.42). These results suggest that hypoxia is positively correlation with m6A/m5C/m1A mediated methylation regulation.

Fig. 1
figure 1

The relationship between hypoxia and m6A/m5C/m1A methylation regulation. A The consensus matrix of unsupervised clustering when K = 2. B Cumulative distribution function (CDF) curve proves that K = 2 has the best clustering effect. C K-M survival curves showed the differences of overall survival rate among the 2 clusters. D Analysis of activation difference of hypoxia pathway between two clusters, the significance of the difference was analyzed by Wilcoxon signed rank test. E Heatmap of differential expression of m6A/m5C/m1A methylation regulation genes in two clusters, the significance of the difference was analyzed by T test. F Analysis of differences in Mean DNA methylation degree between two clusters, the significance of the difference was analyzed by Wilcoxon signed rank test. G Analysis of activation differences of m6A/m5C/m1A methylation regulation between two clusters, the significance of the difference was analyzed by Wilcoxon signed rank test. H Correlation analysis between hypoxia pathway and m6A/m5C/m1A methylation regulation. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, NS means non-significant)

Full size image

Construction of HMRs model of liver cancer

A total of 34 genes related to hypoxia, m6A/m5C/m1A methylation regulation, and differentially expressed in tumors were used as candidate genes for model creation (Fig. 2A). 13 genes related to prognosis were obtained through Univariate Cox regression analysis (Fig. 2B), and candidate genes were screened using LASSO Cox regression (Fig. 2C, D, E) and Multivariate Cox regression (Fig. 2F). Finally, a liver cancer risk score model related to hypoxia and m6A/m5C/m1A methylation regulation was obtained. The risk score was estimated as follows:

Fig. 2
figure 2

Construction of HMRs model of liver cancer A Venn diagram displaying the 34 intersection genes of the hypoxia-associated DEGs and the m6A/m5C/m1A methylation regulation associated DEGs. B Forest plots shows the 13 survival related genes obtained by Univariate Cox regression; P value corrected using B-H correction. C LASSO deviance profile of the 13 genes. D LASSO regression coefficient profile of the 13 genes. E 5 genes coefficient after LASSO analysis. F Final 3 genes forest plot after Multivariate Cox regression

Full size image
$$\mathrm{Risk~score }= (0.152 \times \mathrm{ expr~of~CEP}55) + (0.366 \times \mathrm{expr~of~DPH}2) + (0.313 \times \mathrm{expr~of~SMS})$$

Validation of the prognostic prediction efficacy of HMRs model in liver cancer

Based on the HMRs model, we calculated the risk scores of patients in each dataset and divided all patients into high-risk and low-risk groups based on the median risk scores in the training set. Survival analysis results showed that patients in the high-risk group had a worse prognosis (Fig. 3A, B, C, p < 0.01) and higher mortality (Fig. 3D-I). t-SNE and PCA analysis showed significant clustering of patients in the high- and low-risk groups (Fig. 3J, K, L for t-SNE, Fig. 3M, N, O for PCA). Risk genes were significantly high expression in the high-risk group (Fig. 3P, Q, R). These results indicated that the HMRs model is a robust prediction classifier for liver cancer prognosis.

Fig. 3
figure 3

Validation of the prognostic prediction efficacy of HMRs model A B C K-M curve of OS of liver cancer patients in high- and low-risk groups. D E F Risk curve for the investigated patients. G H I Risk scatter diagram for the investigated patients. The investigated patients risk scores clustering based on t-SNE J K L and PCA M N O P Q R Heatmap for the different expression of 3 risk genes in high- and low-risk group, the significance of the difference was analyzed by T test, the red asterisk represents significantly higher expression in the latter group. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, NS means non-significant)

Full size image

The time-dependent ROC curve revealed that the HMRs model had good predictive performance in both the training and validation sets. The area under the curve (AUC) of 1-year, 2-year, and 3-year of the training set and validation set were greater than 0.65 (0.765, 0.686, 0.689 in TCGA-LIHC, 0.681, 0.711, 0.653 in NODE-HCC, 0.719, 0.677, 0.724 in ICGC-LIHC-JP) (Fig. 4A-C). The results of prognostic independence analysis based on Univariate and Multivariate Cox regression indicated that the HMRs model was an independent prognostic factor for liver cancer (Fig. 4D-I). Univariate Cox regression analysis identified clinical stage and risk score as risk factors for liver cancer (p < 0.001 in all three data sets). Multivariate Cox regression analysis showed that risk score was the only independent prognostic factor in TCGA-LIHC and NODE-HCC datasets. Analysis of the ICGC-LIHC-JP dataset revealed clinical stage, risk score, and gender as risk factors for liver cancer (p < 0.05). These results strongly suggest that HMRs represent an independent predictor of the prognosis of liver cancer patients.

Fig. 4
figure 4

Validation of the prognostic independence of the HMRs model A B C The 1-, 2-, 3-year AUC in ROC analysis. D E F Univariate and G H I Multivariate Cox regression analysis of risk score and clinical features

Full size image

The HMRs model is significantly correlated with genomic alterations

Due to the significant correlation between the HMRs model and the prognosis of liver cancer patients, we conducted a series of bioinformatics analyses to reveal the molecular mechanisms associated with this model.

A total of 893 genes (479 down-regulated, 414 up-regulated) were differentially expressed between high- and low-risk groups (Fig. S2). Pathway analysis showed that DEGs were significantly enriched in processes such as metabolism (monocarboxylic acid metabolic process, DNA metabolic process) and immune response (adaptive immune system) (Fig. 5A). Activation analysis of typical tumor pathways showed that pathways related to cancer progression, such as Androgen, EGFR, Estrogen, Hypoxia, JAK-STAT, and MARK were significantly activated in the high-risk group (Fig. 5B). However, pathways related to cancer suppression and immune response, such as NFKB, p53, TGFβ, TNFα, and WNT, were also significantly activated in the high-risk group (Fig. 5B). Pathways related to cell proliferation and migration, such as PI3K and VEGF, were significantly activated in the low-risk group (Fig. 5B). These results indicate a correlation between the HMRs model and pathways such as cell proliferation and immune regulation.

Fig. 5
figure 5

Relationship between risk models and genomic alterations A Pathway annotation of DEGs between risk groups. B Analysis of differences in activation of 14 typical cancer-related pathways in the high- and low-risk group, the significance of the difference was analyzed by Wilcoxon signed rank test. C Analysis of differences in TMB between two risk groups, the significance of the difference was analyzed by Wilcoxon signed rank test. D E The mutation landscape of patients in the TCGA-LIHC cohort (show the top 20 genes) about the high- and low-risk groups. F Display of mutation genes information with significant differences between risk groups. G Display of mutation sites of TP53 differences between high- and low-risk groups. Display CNV information of high H and low I risk groups. J Analysis of differences in fraction genome altered between two risk groups, the significance of the difference was analyzed by Wilcoxon signed rank test. K Analysis of differences in HRD score between two risk groups, the significance of the difference was analyzed by Wilcoxon signed rank test. Analysis of tumor cell stemness index based on methylation L and epigenetic regulation characteristics M between two risk groups, the significance of the difference was analyzed by Wilcoxon signed rank test. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, NS means non-significant)

Full size image

We then analyzed the differences in gene mutations between high- and low-risk groups. Although there was no significant difference in TMB between the two groups (Fig. 5C), there were differences in mutation landscape between the high- and low-risk groups (Fig. 5D, E). Genes such as TP53, COL7A1, DNAH10, CACNA1H, DIP2C, and NLRP2 exhibited a higher mutation frequency in the high-risk group, while NLRP12 has a higher mutation frequency in the low-risk group (Fig. 5F). Mutation site analysis showed that the type and number of mutations in TP53 were significantly higher than those in the low-risk group (Fig. 5G). These results suggest the presence of differences in gene mutation patterns among patients in the high- and low-risk groups, and there is significant heterogeneity between the two groups of patients. Genomic variation analysis showed that the high-risk group had more CNV and structural variations (Fig. 5H, I, J). The higher HRD score in the high-risk group indicated that there were more defects in the DNA damage repair mechanism in the high-risk group, leading to more mutations (Fig. 5K). The cancer stem cell index based on methylation suggested that the stemness characteristics of the high-risk group were more obvious (Fig. 5L), while the results based on the epigenetic regulation characteristics of stem cells suggested that there was no significant difference between the two groups (Fig. 5M). These results indicate significant differences in expression patterns and gene mutations between high- and low-risk groups, leading to survival differences between groups.

Analysis of immune microenvironment between risk groups

Immune infiltration analysis using the estimate suggested a higher degree of immune infiltration in the high-risk group (Fig. 6A), while there were no significant differences in tumor cell purity, stromal cell score, and estimate score (Fig. S3A, B, C). Immune cell infiltration analyses of 28 types of immune cells showed that Activated CD4 T cells, Activated dendritic cells, Central memory CD4 T cells, Central memory CD8 T cells, Effector memory CD4 T cells, MDSC, Memory B cells, Regulatory T cells, T follicular helper cells, Type 17 T helper cells, Type 2 T helper cells, were significantly activated in the high-risk group, while Eosinophils and CD56dim natural killer cells were highly infiltrated in the low-risk group (Fig. 6B). Other immune infiltration analysis approaches consistently suggested that the high-risk group has more abundant immune cell infiltration (Fig. S4A, B). The infiltrating immune cells in the high-risk group consisted of multiple types of immunosuppressive and activated cells, indicating a more complex immune microenvironment in the high-risk group. TIP analysis suggested that there was no significant difference in immune activation between high- and low-risk groups (Fig. 6C). Analysis of anti-tumor immune pathway activation showed that the high-risk group elevated activation levels in pathway associated with the Release of cancer cell antigens, CD8 T cell recruiting, Neutrophil recruiting, TH17 cell recruiting, and MDSC recruiting pathways in steps 1 and 4, while the Infiltration of immune cells into tumors, Recognition of cancer cells by T cells, Killing of cancer cells in steps 5, 6, 7 were significantly activated in the low-risk group (Fig. 6D). Taken together, these results indicated that the low-risk group could more effectively eliminate tumor cells, thereby improving the prognosis of patients.

Fig. 6
figure 6

Analysis of tumor immune microenvironment A Analysis of differences in immune scores between two risk groups B Analysis of immune cell infiltration differences between two groups. C Overall immune activation differences between two risk groups based on TIP. D Differences in the cancer immunity cycle activities between two groups. The significance of the difference was analyzed by Wilcoxon signed rank test. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, NS means non-significant)

Full size image

Prediction of treatment outcomes for liver cancer using the HMRs model

Analysis of the IMvigor210-based immunotherapy dataset showed that patients who responded to treatment had higher risk scores (Fig. 7A). Evaluation of immunotherapy effectiveness based on TIDE suggested that the high-risk group had a lower TIDE score, indicating a better response to immunotherapy (Fig. 7B). TIS is considered a predictive biomarker for the combination therapy of PD-1 inhibitors, and patients with high TIS scores can benefit from treatment [56, 57]. Analysis showed that the TIS of the high-risk group was higher, indicating that immunotherapy combination therapy conferred more substantial benefits for the high-risk group (Fig. 7C). IPS analysis showed that the high-risk group had significantly higher IPS scores (Fig. 7D), while the remaining scores were comparable to the low-risk group (Fig. 7E, F, G). These results indicate that the high-risk group was more likely to benefit from immunotherapy. Multiple common ICIs and HLA were significantly expressed in the high-risk group (Fig. 7H, I), which also supports the speculation that the high-risk group can derive more benefit from immunotherapy. We also analyzed the predictive ability of our risk score for the effectiveness of TACE treatment. The results showed that TACE had a better effect on the low-risk group (Fig. 7J). Finally, we analyzed the guidance of risk scoring on chemotherapy medication. Of these 8 commonly used chemotherapy drugs, Erlotinib was more suitable for the low-risk group, while Cisplatin, 5-Fluorouracil, Vorinostat, Tivozanib, Doxorubicin, and Temsirolimus were more beneficial to the high-risk group (Fig. 7K).

Fig. 7
figure 7

The ability of HMRs model to predict the outcome of liver cancer treatment A Association between risk scores and immunotherapy response based on the IMvigor210. B Difference in response to immunotherapy in risk groups predicted by TIDE. Differences in the TIS C IPS D IPS-CTLA4- and PD1/PDL1/PDL2 blocker E IPS-CTLA4 blocker F, IPS-PD1/PDL1/PDL2 blocker G ICIs H and HLA I between two groups. J Association between risk scores and TACE treatment response. K Difference in IC50 of common chemotherapy drugs between risk group. The significance of the difference was analyzed by Wilcoxon signed rank test. (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, NS means non-significant)

Full size image

Analysis of expression patterns of risk genes

Due to the lack of clinical samples for risk gene validation in this study, we analyzed the expression patterns of risk genes from multiple perspectives. Analysis based on the ICGC-JP, GSE144269, and TCGA-LIHC datasets showed significant overexpression of risk genes in tumor tissues (Fig. S5A, B, C). Analysis based on the UALCAN database showed that the proteins of DPH2 and SMS were significantly overexpressed in tumors, while the difference in CEP55 was not significant (Fig. S5D). The results of immunohistochemistry analysis showed that CEP55 was moderately expressed in tumor tissue, with a lower expression level in normal tissue (Fig. S5E), while SMS exhibited low expression in tumor tissue and could not be detected in normal tissue (Fig. S5G). DPH2 was highly expressed in tumor tissues, while its expression level was low in normal tissues (Fig. S5F). We validated the expression patterns of risk genes using the expression chip data collected by TNM plot for liver cancer. The results showed that risk genes were significantly overexpressed in tumor tissues, both in paired and non-paired samples (Fig. S5H, I, J). In addition, pan-cancer analysis based on GEPIA2 showed that high expression of risk genes was associated with poor prognosis in various cancers (Fig. S6A, B, C, D, including liver cancer). Based on the above results, risk genes are significantly overexpressed in liver cancer’s tumor tissue and are all risk factors for liver cancer.

Gene function analysis of DPH2

The results of multivariate Cox regression analysis highlighted the significant impact of DPH2 on the prognosis of liver cancer (p = 0.014, hazard ratio > 1, Fig. 2F). We also confirmed that DPH2 was significantly overexpressed in both gene and protein expression levels in liver tumor tissue. Therefore, we analyzed the role of this gene in the development of liver cancer. Univariate and Multivariate Cox regression analysis showed that DPH2 was a significant independent prognostic factor in both the TCGA-LIHC (Fig. S7A, D) and ICGC-LIHC-JP (Fig. S7C, F) datasets, but nor significant in the NODE-HCC dataset (Fig. S7B, E). The time-dependent ROC curve analysis of DPH2 expression suggested that the AUC values of both TCGA-LIHC (Fig. S7G) and ICGC-LIHC-JP (Fig. S7I) datasets for survival at 1, 2, and 3 years are greater than 0.65, indicating that the gene has good predictive ability for liver cancer prognosis, while the NODE-HCC dataset revealed that DPH2 lacks predictive ability (Fig. S7H). Analysis shows that normal and low-level tumor-grade samples have higher levels of DPH2 methylation and lower gene expression (Fig. S7J, L). Patients with low gene expression levels were associated with higher methylation levels (Fig. S7K). Multiple methylation modification sites were observed in the promoter region of the DPH2 gene (Fig. S7M). These results suggest that DPH2 is a potential prognostic target in liver cancer, and methylation in the promoter region of this gene may affect its expression, thereby influencing clinical outcomes.

Further analysis showed that the correlation between DPH2 and immune cells is not significant (Fig. S6E). Analysis of liver cancer single-cell data based on TISCH2 showed that DPH2 exhibits low expression in immune-related cells (Fig. S8A). PanglaoDB database analysis showed that DPH2 is highly expressed in Basal cells (Fig. S8B). THPA database analysis demonstrated that the protein expression level of DPH2 increased with the progression of cell interphase, unlike gene expression (Fig. S8C, D, E, F). These results suggest a possible regulatory relationship between DPH2 and cell proliferation.

Discussion

Liver cancer exhibits significant inter-tumor heterogeneous and involves complex and diverse biological processes, resulting in inconsistent prognoses [3, 68]. Given its high morbidity and mortality rates and limited treatment options, it is imperative to develop a convenient method for guiding personalized precision treatment for liver cancer.

Current evidence suggests that hypoxia within the TME is closely associated with unfavorable processes in liver cancer, including tumor proliferation, metastasis, angiogenesis, resistance to radio- and chemotherapy, and disease progression [69]. The impact of hypoxia on the TME is multifaceted, encompassing metabolic reprogramming and the modulation of immune cell functions through post-transcriptional methylation modifications [28]. Notably, there exists an intricate crosstalk between hypoxia in solid tumors and methylation modifications such as m6A/m5C/m1A [15]. In this study, we comprehensively examined the relationship between hypoxia and m6A/m5C/m1A methylation and subsequently developed a prognostic risk model for liver cancer by integrating hypoxia and methylation information. To our knowledge, this represents the first reported liver cancer prognosis model that integrates gene expression related to hypoxia and methylation regulation. Based on our comprehensive analyses, the HMRs model serves as a robust classifier for predicting liver cancer prognosis and offers valuable insights for guiding liver cancer immunotherapy and chemotherapy.

Our analysis preliminarily revealed a positive correlation between hypoxia and methylation in liver cancer, as patients with higher hypoxia pathway activation exhibited poorer prognoses, along with elevated methylation levels and increase m6A/m5C/m1A pathway activation. These finding validate the rationale behind our integrated modeling approach. Survival curves, clustering, ROC analysis, and both Univariate and Multivariate Cox analyses of the training and validation sets provided compelling evidence of the effectiveness and robustness of HMRs in predicting liver cancer prognosis.

Hypoxia has been linked to the promotion of tumor heterogeneity and genomic evolution, with high-hypoxia patients often displaying elevated CNV and TMB [70]. Additionally, hypoxia contributes to reshaping the TME by inducing m6A/m5C/m1A methylation, influencing biological functions such as immunosuppression, metabolic dysregulation, and metastasis promotion [15]. Studies have shown that genomic instability represented by TMB and abnormal immune microenvironment can reflect tumor prognosis and immunotherapy effect [71, 72]. In our study, we conducted a series of analyses using the TCGA-LIHC dataset to investigate the biological processes associated with HMRs. Our analyses revealed significant disparities in immune responses, cell proliferation, and metabolism pathways between high- and low-risk groups. Mutation of TP53 will affect cell apoptosis and cause the tumor to fall into a dynamic cycle of hypoxia [73]. Although no significant differences in TMB were observed between risk groups, the mutated genes were notably different, with TP53 mutations being more prevalent in the high-risk group, which also exhibited a higher frequency of CNV. Furthermore, we observed increased HRD in the high-risk group, suggesting greater genomic instability. Notably, hypoxia’s impact on gene expression, especially that of m6A/m5C/m1A regulators, has been shown to enrich and maintain tumor stem cells, consistent with our observation of a higher cancer stem cell index in the high-risk group [74, 75]. These finding collectively indicate that the HMRs model encompasses the characteristics of hypoxia and m6A/m5C/m1A methylation in regulating cancer biological processes, making it a potential prognostic predictor for liver cancer.

Further analysis revealed that the low-risk group exhibit a higher proportion of Eosinophils and CD56dim natural killer cells, potentially contributing to cancer cell suppression and improved survival. Eosinophils are thought to enhance immune function, while CD56dim natural killer cells inhibit cancer cell proliferation [76, 77]. Moreover, our assessment of anti-tumor immune pathway activation using the TIP tool suggested that the low-risk group may eliminate tumor cells more efficiently, contributing to better patient prognoses. Conversely, the high-risk group display complex immune cell infiltration with both immune-activating and immunosuppressive cell types exhibiting higher proportions. While immune activating cells like Central memory CD4 T cells, Central memory CD8 T cells and Type 2 T helper cells were more abundant in the high-risk group [78,79,80]. Immunosuppressive or inflammation-promoting cells like Effector memory CD4 T cells, MDSC and Regulatory T cells were also highly prevalent in the high-risk group [81, 82]. The effect of immune cells in TME on the occurrence and metastasis of hepatocellular carcinoma is complex [83]. This intricate immune cell composition in the high-risk group may renders anti-tumor pathways less effective, resulting in poorer prognoses.

Our results also indicated that the HMRs model could effectively assess the efficacy of immunotherapy and TACE for liver cancer. Despite limited immunotherapy datasets for liver cancer, multiple analysis methods consistently demonstrated that the high-risk group display better responses to immunotherapy. The high group’s increased release of cancer cell antigens and higher expression of multiple ICIs and HLAs suggest a potential advantage in responding to immunotherapy, supported by the release of more tumor antigens. Studies have shown that combined radiotherapy, chemotherapeutic, and immunotherapy can improve the prognosis of patients with liver cancer [84]. Our study also identified specific chemotherapeutic drugs that showed varying responses between risk groups, offering potential clinical benefits, particularly when combined with immunotherapy. TACE, on the other hand, exhibited better responses in the low -risk group, thereby providing valuable treatment guidance for patients in this subgroup.

Furthermore, we investigated the gene DPH2, known to be involved in ribosomal protein biosynthesis and protein elongation accuracy [85]. Expression analyses across multiple datasets consistently indicated significant overexpression of DPH2 in liver cancer. We observed a positive correlation between DPH2 expression levels and tumor stage, coupled with a negative correlation between methylation levels and both DPH2 expression and clinical phenotypes. These findings suggest a close association between DPH2 and mitosis. In normal tissues, promoter hypermethylation keeps DPH2 expression low, while in cancer tissues, decreasing methylation levels with increasing tumor stage elevate DPH2 expression potentially influencing disease progression. We propose that DPH2 represents a potential target for liver cancer prognosis evaluation and treatment.

Conclusion

In summary, our study unveiled the intricate relationship between hypoxia and m6A/m5C/m1A methylation regulation, leading to the development of the Hypoxia-m6A/m5C/m1A-related-Score model. This model captures the biological phenotypes of hypoxia and m6A/m5C/m1A methylation, such as genomic instability, metabolic abnormalities, and immune dysregulation, effectively predicting liver cancer survival outcomes. Moreover, our model offers valuable insight into guiding immunotherapy, TACE, and chemotherapy for liver cancer, thus facilitating precise treatment strategies for liver cancer patients.

Availability of data and materials

The original contributions presented in the study are included in the article/Supplement Material.

References

  1. Cisneros-Garza LE, González-Huezo MS, Moctezuma-Velázquez C, Ladrón de Guevara-Cetina L, Vilatobá M, García-Juárez I, Alvarado-Reyes R, Álvarez-Treviño GA, Allende-Pérez S, Bornstein-Quevedo L, et al. The second Mexican consensus on hepatocellular carcinoma. Part I: Epidemiology and diagnosis. Rev Gastroenterol Mex (Engl Ed). 2022;87:216–34.

    CAS  PubMed  Google Scholar 

  2. Yang JD, Heimbach JK. New advances in the diagnosis and management of hepatocellular carcinoma. BMJ. 2020;371:m3544.

    Article  PubMed  Google Scholar 

  3. Craig AJ, von Felden J, Garcia-Lezana T, Sarcognato S, Villanueva A. Tumour evolution in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. 2020;17:139–52.

    Article  PubMed  Google Scholar 

  4. Li H, Li CW, Li X, Ding Q, Guo L, Liu S, Liu C, Lai CC, Hsu JM, Dong Q, et al. MET inhibitors promote liver tumor evasion of the immune response by stabilizing PDL1. Gastroenterology. 2019;156:1849–1861.e1813.

    Article  CAS  PubMed  Google Scholar 

  5. Al-Salama ZT, Syed YY, Scott LJ. Lenvatinib: a review in hepatocellular carcinoma. Drugs. 2019;79:665–74.

    Article  CAS  PubMed  Google Scholar 

  6. Graham K, Unger E. Overcoming tumor hypoxia as a barrier to radiotherapy, chemotherapy and immunotherapy in cancer treatment. Int J Nanomedicine. 2018;13:6049–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Daşu A, Toma-Daşu I, Karlsson M. Theoretical simulation of tumour oxygenation and results from acute and chronic hypoxia. Phys Med Biol. 2003;48:2829–42.

    Article  PubMed  Google Scholar 

  8. Gilkes DM, Semenza GL, Wirtz D. Hypoxia and the extracellular matrix: drivers of tumour metastasis. Nat Rev Cancer. 2014;14:430–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Bertout JA, Patel SA, Simon MC. The impact of O2 availability on human cancer. Nat Rev Cancer. 2008;8:967–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Parks SK, Cormerais Y, Pouyssegur J. Hypoxia and cellular metabolism in tumour pathophysiology. J Physiol. 2017;595:2439–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Horsman MR, Mortensen LS, Petersen JB, Busk M, Overgaard J. Imaging hypoxia to improve radiotherapy outcome. Nat Rev Clin Oncol. 2012;9:674–87.

    Article  CAS  PubMed  Google Scholar 

  12. Jing X, Yang F, Shao C, Wei K, Xie M, Shen H, Shu Y. Role of hypoxia in cancer therapy by regulating the tumor microenvironment. Mol Cancer. 2019;18:157.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Shannon AM, Bouchier-Hayes DJ, Condron CM, Toomey D. Tumour hypoxia, chemotherapeutic resistance and hypoxia-related therapies. Cancer Treat Rev. 2003;29:297–307.

    Article  CAS  PubMed  Google Scholar 

  14. Schlageter M, Terracciano LM, D’Angelo S, Sorrentino P. Histopathology of hepatocellular carcinoma. World J Gastroenterol. 2014;20:15955–64.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Zhang F, Liu H, Duan M, Wang G, Zhang Z, Wang Y, Qian Y, Yang Z, Jiang X. Crosstalk among m(6)A RNA methylation, hypoxia and metabolic reprogramming in TME: from immunosuppressive microenvironment to clinical application. J Hematol Oncol. 2022;15:84.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Teng PC, Liang Y, Yarmishyn AA, Hsiao YJ, Lin TY, Lin TW, Teng YC, Yang YP, Wang ML, Chien CS, et al. RNA modifications and epigenetics in modulation of lung cancer and pulmonary diseases. Int J Mol Sci. 2021;22(19):10592.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Yang Z, Wang T, Wu D, Min Z, Tan J, Yu B. RNA N6-methyladenosine reader IGF2BP3 regulates cell cycle and angiogenesis in colon cancer. J Exp Clin Cancer Res. 2020;39:203.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Zhao Y, Zhao Q, Kaboli PJ, Shen J, Li M, Wu X, Yin J, Zhang H, Wu Y, Lin L, et al. m1A regulated genes modulate PI3K/AKT/mTOR and ErbB pathways in gastrointestinal cancer. Transl Oncol. 2019;12:1323–33.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Yang H, Wang Y, Xiang Y, Yadav T, Ouyang J, Phoon L, Zhu X, Shi Y, Zou L, Lan L. FMRP promotes transcription-coupled homologous recombination via facilitating TET1-mediated m5C RNA modification demethylation. Proc Natl Acad Sci U S A. 2022;119: e2116251119.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Harris AL. Hypoxia–a key regulatory factor in tumour growth. Nat Rev Cancer. 2002;2:38–47.

    Article  CAS  PubMed  Google Scholar 

  21. Wigerup C, Påhlman S, Bexell D. Therapeutic targeting of hypoxia and hypoxia-inducible factors in cancer. Pharmacol Ther. 2016;164:152–69.

    Article  CAS  PubMed  Google Scholar 

  22. Hu CJ, Wang LY, Chodosh LA, Keith B, Simon MC. Differential roles of hypoxia-inducible factor 1alpha (HIF-1alpha) and HIF-2alpha in hypoxic gene regulation. Mol Cell Biol. 2003;23:9361–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Liu Y, Yan W, Tohme S, Chen M, Fu Y, Tian D, Lotze M, Tang D, Tsung A. Hypoxia induced HMGB1 and mitochondrial DNA interactions mediate tumor growth in hepatocellular carcinoma through Toll-like receptor 9. J Hepatol. 2015;63:114–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Daniel SK, Sullivan KM, Labadie KP, Pillarisetty VG. Hypoxia as a barrier to immunotherapy in pancreatic adenocarcinoma. Clin Transl Med. 2019;8:10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Chouaib S, Noman MZ, Kosmatopoulos K, Curran MA. Hypoxic stress: obstacles and opportunities for innovative immunotherapy of cancer. Oncogene. 2017;36:439–45.

    Article  CAS  PubMed  Google Scholar 

  26. Dong F, Qin X, Wang B, Li Q, Hu J, Cheng X, Guo D, Cheng F, Fang C, Tan Y, et al. ALKBH5 facilitates hypoxia-induced paraspeckle assembly and IL8 Secretion to generate an immunosuppressive tumor microenvironment. Cancer Res. 2021;81:5876–88.

    Article  CAS  PubMed  Google Scholar 

  27. Yang N, Wang T, Li Q, Han F, Wang Z, Zhu R, Zhou J. HBXIP drives metabolic reprogramming in hepatocellular carcinoma cells via METTL3-mediated m6A modification of HIF-1α. J Cell Physiol. 2021;236:3863–80.

    Article  CAS  PubMed  Google Scholar 

  28. Ruan DY, Li T, Wang YN, Meng Q, Li Y, Yu K, Wang M, Lin JF, Luo LZ, Wang DS, et al. FTO downregulation mediated by hypoxia facilitates colorectal cancer metastasis. Oncogene. 2021;40:5168–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Damgaci S, Ibrahim-Hashim A, Enriquez-Navas PM, Pilon-Thomas S, Guvenis A, Gillies RJ. Hypoxia and acidosis: immune suppressors and therapeutic targets. Immunology. 2018;154:354–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Wang B, Zhao Q, Zhang Y, Liu Z, Zheng Z, Liu S, Meng L, Xin Y, Jiang X. Targeting hypoxia in the tumor microenvironment: a potential strategy to improve cancer immunotherapy. J Exp Clin Cancer Res. 2021;40:24.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics. 2007;23:1846–7.

    Article  PubMed  Google Scholar 

  32. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Wang E, Li Y, Ming R, Wei J, Du P, Zhou P, Zong S, Xiao H. The Prognostic Value and Immune Landscapes of a m(6)A/m(5)C/m(1)A-Related LncRNAs Signature in Head and Neck Squamous Cell Carcinoma. Front Cell Dev Biol. 2021;9:718974.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Li D, Li K, Zhang W, Yang KW, Mu DA, Jiang GJ, Shi RS, Ke D. The m6A/m5C/m1A regulated gene signature predicts the prognosis and correlates with the immune status of hepatocellular carcinoma. Front Immunol. 2022;13:918140.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26:1572–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10:1523.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28:1747–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12:R41.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Ramos M, Geistlinger L, Oh S, Schiffer L, Azhar R, Kodali H, de Bruijn I, Gao J, Carey VJ, Morgan M, Waldron L. Multiomic integration of public oncology databases in bioconductor. JCO Clin Cancer Inform. 2020;4:958–71.

    Article  PubMed  Google Scholar 

  41. Knijnenburg TA, Wang L, Zimmermann MT, Chambwe N, Gao GF, Cherniack AD, Fan H, Shen H, Way GP, Greene CS, et al. Genomic and molecular landscape of DNA damage repair deficiency across the cancer genome atlas. Cell Rep. 2018;23:239–254.e236.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Walcher L, Kistenmacher AK, Suo H, Kitte R, Dluczek S, Strauß A, Blaudszun AR, Yevsa T, Fricke S, Kossatz-Boehlert U. Cancer stem cells-origins and biomarkers: perspectives for targeted personalized therapies. Front Immunol. 2020;11:1280.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, Kamińska B, Huelsken J, Omberg L, Gevaert O, et al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell. 2018;173:338–354.e315.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. de Visser KE, Joyce JA. The evolving tumor microenvironment: from cancer initiation to metastatic outgrowth. Cancer Cell. 2023;41:374–403.

    Article  PubMed  Google Scholar 

  45. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    Article  PubMed  Google Scholar 

  46. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48:W509–w514.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. 2018;1711:243–59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 2017;18:220.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Plattner C, Finotello F, Rieder D. Deconvoluting tumor-infiltrating immune cells from RNA-seq data using quanTIseq. Methods Enzymol. 2020;636:261–85.

    Article  CAS  PubMed  Google Scholar 

  50. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, Selves J, Laurent-Puig P, Sautès-Fridman C, Fridman WH, de Reyniès A. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016;17:218.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Racle J, Jonge K, Baumgaertner P, Speiser DE, Gfeller D. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife. 2017;6:e26476.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Sturm G, Finotello F, Petitprez F, Zhang JD, Baumbach J, Fridman WH, List M, Aneichyk T. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics. 2019;35:i436–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, Yuan H, Cheng P, Li F, Long Z, et al. TIP: a web server for resolving tumor immunophenotype profiling. Cancer Res. 2018;78:6575–80.

    Article  CAS  PubMed  Google Scholar 

  54. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554:544–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Fu J, Li K, Zhang W, Wan C, Zhang J, Jiang P, Liu XS. Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 2020;12:21.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Danaher P, Warren S, Lu R, Samayoa J, Sullivan A, Pekker I, Wallden B, Marincola FM, Cesano A. Pan-cancer adaptive immune resistance as defined by the Tumor Inflammation Signature (TIS): results from The Cancer Genome Atlas (TCGA). J Immunother Cancer. 2018;6:63.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Ayers M, Lunceford J, Nebozhyn M, Murphy E, Loboda A, Kaufman DR, Albright A, Cheng JD, Kang SP, Shankaran V, et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. 2017;127:2930–40.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 2017;18:248–62.

    Article  CAS  PubMed  Google Scholar 

  59. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS ONE. 2014;9:e107468.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, Bindal N, Beare D, Smith JA, Thompson IR, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 2013;41:D955–961.

    Article  CAS  PubMed  Google Scholar 

  61. Chandrashekar DS, Karthikeyan SK, Korla PK, Patel H, Shovon AR, Athar M, Netto GJ, Qin ZS, Kumar S, Manne U, et al. UALCAN: an update to the integrated cancer data analysis platform. Neoplasia. 2022;25:18–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Sjöstedt E, Zhong W, Fagerberg L, Karlsson M, Mitsios N, Adori C, Oksvold P, Edfors F, Limiszewska A, Hikmet F, et al. An atlas of the protein-coding genes in the human, pig, and mouse brain. Science. 2020;367(6482):eaay5947.

    Article  PubMed  Google Scholar 

  63. Bartha Á, Győrffy B. TNMplot.com: A web tool for the comparison of gene expression in normal, tumor and metastatic tissues. Int J Mol Sci. 2021;22(5):2622.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Tang Z, Kang B, Li C, Chen T, Zhang Z. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 2019;47:W556–w560.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Ou J, Zhu LJ. trackViewer: a Bioconductor package for interactive and integrative visualization of multi-omics data. Nat Methods. 2019;16:453–4.

    Article  CAS  PubMed  Google Scholar 

  66. Sun D, Wang J, Han Y, Dong X, Ge J, Zheng R, Shi X, Wang B, Li Z, Ren P, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. 2021;49:D1420–d1430.

    Article  CAS  PubMed  Google Scholar 

  67. Franzén O, Gan LM, Björkegren JLM. PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data. Database (Oxford). 2019;2019:baz046.

    Article  PubMed  Google Scholar 

  68. Ho DW, Tsui YM, Sze KM, Chan LK, Cheung TT, Lee E, Sham PC, Tsui SK, Lee TK, Ng IO. Single-cell transcriptomics reveals the landscape of intra-tumoral heterogeneity and stemness-related subpopulations in liver cancer. Cancer Lett. 2019;459:176–85.

    Article  CAS  PubMed  Google Scholar 

  69. Wu XZ, Xie GR, Chen D. Hypoxia and hepatocellular carcinoma: The therapeutic target for hepatocellular carcinoma. J Gastroenterol Hepatol. 2007;22:1178–82.

    Article  CAS  PubMed  Google Scholar 

  70. Bhandari V, Li CH, Bristow RG, Boutros PC. Divergent mutational processes distinguish hypoxic and normoxic tumours. Nat Commun. 2020;11:737.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Xue W, Wang Y, Xie Y, Yang C, Gong Z, Guan C, Wei C, Zhu C, Niu Z. miRNA-based signature associated with tumor mutational burden in colon adenocarcinoma. Front Oncol. 2021;11:634841.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Xue W, Dong B, Wang Y, Xie Y, Li P, Gong Z, Niu Z. A novel prognostic index of stomach adenocarcinoma based on immunogenomic landscape analysis and immunotherapy options. Exp Mol Pathol. 2022;128:104832.

    Article  CAS  PubMed  Google Scholar 

  73. Greijer AE, van der Wall E. The role of hypoxia inducible factor 1 (HIF-1) in hypoxia induced apoptosis. J Clin Pathol. 2004;57:1009–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Wang X, Ji Y, Feng P, Liu R, Li G, Zheng J, Xue Y, Wei Y, Ji C, Chen D, Li J. The m6A reader IGF2BP2 regulates macrophage phenotypic activation and inflammatory diseases by stabilizing TSC1 and PPARγ. Adv Sci (Weinh). 2021;8:2100209.

    Article  CAS  PubMed  Google Scholar 

  75. Zhang C, Samanta D, Lu H, Bullen JW, Zhang H, Chen I, He X, Semenza GL. Hypoxia induces the breast cancer stem cell phenotype by HIF-dependent and ALKBH5-mediated m6A-demethylation of NANOG mRNA. Proc Natl Acad Sci U S A. 2016;113:E2047–2056.

    CAS  PubMed  PubMed Central  Google Scholar 

  76. Rosenberg HF, Dyer KD, Foster PS. Eosinophils: changing perspectives in health and disease. Nat Rev Immunol. 2013;13:9–22.

    Article  CAS  PubMed  Google Scholar 

  77. Crinier A, Narni-Mancinelli E, Ugolini S, Vivier E. SnapShot: natural killer cells. Cell. 2020;180:1280–1280.e1281.

    Article  CAS  PubMed  Google Scholar 

  78. Künzli M, Masopust D. CD4(+) T cell memory. Nat Immunol. 2023;24:903–14.

    Article  PubMed  PubMed Central  Google Scholar 

  79. Turner SJ, Bennett TJ, La Gruta NL. CD8(+) T-Cell Memory: The Why, the When, and the How. Cold Spring Harb Perspect Biol. 2021;13(5):a038661.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Walker JA, McKenzie ANJ. T(H)2 cell development and function. Nat Rev Immunol. 2018;18:121–33.

    Article  CAS  PubMed  Google Scholar 

  81. McDaniel MM, Chawla AS, Jain A, Meibers HE, Saha I, Gao Y, Jain V, Roskin K, Way SS, Pasare C. Effector memory CD4(+) T cells induce damaging innate inflammation and autoimmune pathology by engaging CD40 and TNFR on myeloid cells. Sci Immunol. 2022;7:eabk0182.

    CAS  PubMed  PubMed Central  Google Scholar 

  82. Hegde S, Leader AM, Merad M. MDSC: Markers, development, states, and unaddressed complexity. Immunity. 2021;54:875–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Chen Y, Zhou Y, Yan Z, Tong P, Xia Q, He K. Effect of infiltrating immune cells in tumor microenvironment on metastasis of hepatocellular carcinoma. Cell Oncol (Dordr). 2023;46(6):1595–604.

    Article  CAS  PubMed  Google Scholar 

  84. Choi C, Yoo GS, Cho WK, Park HC. Optimizing radiotherapy with immune checkpoint blockade in hepatocellular carcinoma. World J Gastroenterol. 2019;25:2416–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Hawer H, Mendelsohn BA, Mayer K, Kung A, Malhotra A, Tuupanen S, Schleit J, Brinkmann U, Schaffrath R. Diphthamide-deficiency syndrome: a novel human developmental disorder and ribosomopathy. Eur J Hum Genet. 2020;28:1497–508.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by the Natural Science Basic Research Project of Shaanxi Province (2022JQ-159).

Author information

Authors and Affiliations

  1. College of Life Science, Yan’an University, 716000, Yan’an, China

    Min ren, Bei Fan, Rongrong Zong, Liaoliao Feng & Huiru Sun

  2. The First Clinical Medical College, Yan’an University, 716000, Yan’an, China

    Guangcai Cao

Authors
  1. Min ren
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Bei Fan
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Guangcai Cao
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Rongrong Zong
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Liaoliao Feng
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Huiru Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

MR and HS designed study, wrote the manuscript, and collected the corresponding datasets. BF, GC and RZ assisted in the model construction and model validation. MR and LF drew the figures as well as tables. HS supervised the whole project and edited the manuscript. All authors read and approve the final manuscript.

Corresponding author

Correspondence to Huiru Sun.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

ren, M., Fan, B., Cao, G. et al. Exploration and validation of a combined Hypoxia and m6A/m5C/m1A regulated gene signature for prognosis prediction of liver cancer. BMC Genomics 24, 776 (2023). https://doi.org/10.1186/s12864-023-09876-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09876-3

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us