Multi-omics landscape of m6A methylation regulators in HCC
The roles of 26 m6A RNA methylation regulating genes in HCC were investigated in this study, including 9 writers (METTL3, METTL14, METTL16 [Also known as METT10D], WTAP, KIAA1429, ZC3H13, RBM15, RBM15B, and CBLL1), 3 erasers (ALKBH3, ALKBH5, and FTO), and 14 readers (YTHDF1, YTHDF2, YTHDF3, YTHDC1, YTHDC2, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, HNRNPC, FMR1, LRPPRC, ELAVL1, and EIF3A). Figure 1A depicts the 26 m6A methylation regulatory factors that can recognize, delete, and add m6A-modified sites through a potential regulatory mechanism. However, it is not clear how m6A methylation factors are specifically and dynamically deposited in RNA.
In the HCC genomes of the TCGA and ICGC cohorts, the total mutation rate of 26 m6A regulators is relatively low. Forty-four samples in the TCGA cohort and 37 samples in the ICGC cohort experienced genetic alterations of m6A regulators (Fig. 1B and Figure S2A, respectively), and the most common alteration was missense mutation. In the TCGA cohort, KIAA1429 had the greatest mutation frequency, followed by HNRNPC and LRPPRC (Fig. 1B), whereas in the ICGC cohort, HNRNPA2B1 had the highest mutation frequency, followed by KIAA1429 and EIF3A (Figure S2A). We also conducted a co-mutation analysis on 26 m6A regulators and discovered co-occurrence mutation patterns in YTHDC2 and EIF3A, IGF2BP2 and KIAA1429, FMR1 and HNRNPC (Figure S2B). However, no significant co-exclusive mutation event was found among 26 m6A regulators (Figure S2B). Further analysis copy number variation (CNV) of 26 m6A regulators revealed that KIAA1429, YTHDF3, and HNRNPA2B1 had a relatively high frequency of amplification, while METTL16, METTL14, and ZC3H13 had frequent CNV deletions (Fig. 1D). Figure S2C illustrates the chromosomal locations of CNV mutations in 26 m6A regulators. Furthermore, we utilize STRING database to evaluate PPI network and found that 26 m6A regulators had widespread protein interactions (Fig. 1C). Most regulators were strongly elevated in HCC tissues, with the exception of ZC3H13 and YTHDC2, as shown in Fig. 1E, implying that they may play essential roles in the formation and progression of HCC. Mutual regulation in the TCGA cohort was assessed via Spearman correlation test. We found that most of these m6A regulators had a strong positive co-expression relationship (Figure S2D), in accord with the results in the gathered HCC cohort (Fig. 2A and Table S4). In conclusion, these findings revealed that genetic variation and abnormal expression in the 26 m6A regulators may play critical roles in HCC initiation and progression.
Prognosis And Immune Characteristics Of Ma Regulators
Four datasets with available survival information and clinical data (TCGA-LIHC, ICGC-LIRI-JP, GSE76427, and GSE116174) were enrolled in our study. The prognostic value of 26 m6A regulators was assessed utilizing univariate Cox regression analysis. Most regulators, such as RBM15B, YTHDF2, and LRPPRC, were recognized as risk factors in a forest plot, whereas ZC3H13 could be regarded as a protective factor, as it was shown to be significantly associated with longer overall survival (Fig. 2A and S3A). The network of 26 m6A regulators provided a comprehensive profile of these regulators' interconnections and prognostic value (Fig. 2A). The findings suggested that cross-talk between these regulators could imply a functional relationship and play a pivotal role in the development and progression of HCC. We explored the interaction between 26 m6A regulators and immune infiltrating cells to better understand the role of these m6A regulators in HCC. WTAP and ALKBH3 expression was positively correlated with most immune cells, whereas the remaining m6A regulators were negatively correlated with most immune cells (Fig. 2B). Then, we performed differential analysis of 26 m6A regulators among the four immune response phenotypes (complete response, CR; partial response, PR; stable disease, SD; progressive disease, PD) in the IMvigor210 cohort (Figure S4). As showed in Fig. 2C, there was a substantial difference in expression of the reader HNRNPA2B1 among the four immune response phenotypes, with CR/PR showing higher expression than SD/PD. The Kaplan-Meier method was used to evaluate the prognostic value of HNRNPA2B1 in HCC. Survival analysis revealed that patients with HNRNPA2B1 high expression had a poor outcome in the TCGA LIHC cohort (Fig. 2D). Conversely, patients with high HNRNPA2B1 expression had a better prognosis in the IMvigor210 immunotherapy cohort (Fig. 2E). The molecular mechanisms underlying the abnormal phenomenon may suggest to a functional correlation that needs to be further revealed. GSEA analysis in the TCGA cohort revealed that several cancer-related and immune-related pathways (such as the FC γ R-mediated phagocytosis and T cell receptor signaling pathway) were markedly enriched in the HNRNPA2B1 high-expression group (Fig. 2F). The findings suggest a preliminary correlation between immunotherapy response, cancer formation and progression, and HNRNPA2B1 expression. The reader HNRNPA2B1 might be a potential biomarker for predicting immunotherapy response.
Identification of m 6 A methylation modification patterns mediated by 26 m 6 A regulators in the gathered HCC cohort
To further understand m6A methylation modification patterns, we deployed the consensus clustering approach to stratify samples in the gathered HCC cohort based on the expression of 26 m6A regulators. As shown in Fig. 3A, we identified two distinct modification clusters that could provide the greatest clustering efficacy. Patients were then divided into two categories: Cluster 1 (n = 461) and Cluster 2 (n = 318). We termed these patterns C1 and C2. Most m6A regulators showed relatively high expression in pattern C2 (Figure S3B, and S3C). Based on the expression of 26 m6A methylation regulators, principal component analysis (PCA) revealed that two subgroups were well separated, suggesting that they were effectively distinguishable (Fig. 3B). In the gathered HCC cohort, survival analysis showed that pattern C1 had a superior survival advantage, whereas pattern C2 had the poor prognosis (P < 0.0001, log-rank test, Fig. 3C). We subsequently explored the relative abundances of 28 immune infiltrating cells between different m6A modification patterns. Table S5 lists the relative abundance of 28 immune cells in TME cell infiltration between patterns C1 and C2. As shown in Fig. 3D, pattern C2 did not exhibit a survival advantage over pattern C1, even though pattern C2 had significantly more activated CD4+ T cells, effector memory CD4+ T cells, and natural killer T cells, as well as fewer regulatory T cells, type 1 T helper cells, and monocytes. Pattern C1 represented favorable overall survival and characterized by a relatively high proportion of activated CD8+ T cells and effector memory CD8+ T cells.
Then, the overall infiltration of immune cells (Immune Score), stromal cells (Stromal Score), and tumor cell purity (Tumor Purity) between the two m6A modification patterns were calculated utilizing ESTIMATE algorithm. Pattern C1 had higher stromal scores than pattern C2 (Fig. 4A, P = 0.00021). Pattern C2 had a greater tumor purity than pattern C1, implying that the pattern C2 subtype contains more tumor components (Fig. 4B, P = 0.0079). The immune score between pattern C1 and C2 did not achieve statistically significant (Fig. 4C, P = 0.096). In addition, the differences between the two m6A methylation regulator patterns in signaling pathways were investigated. The m6A methylation pattern C1 was significantly enriched in metabolism-related pathways, such as xenobiotic metabolism, bile acid metabolism, and fatty acid metabolism, etc., whereas pattern C2 was markedly enhanced in cell cycle-related pathways, as exhibited by processes related to MYC targets V1/V2, DNA repair, and E2F target pathways, indicating a cancer activation status in pattern C2 (Fig. 4D, and Table S6).
Different Clinical And Transcriptome Characteristics Of The Two Ma Methylation Patterns In The Tcga-lihc Cohort
The TCGA-LIHC cohort was employed to further validate putative m6A methylation patterns in HCC. Consensus clustering of the expression of 26 m6A regulators also revealed two m6A modification patterns (Figure S5A). The two patterns of 26 m6A methylation regulators in the TCGA LIHC cohort were almost similar to that of the gathered HCC cohort. In accord with the results in the gathered HCC cohort, a similar expression pattern of 26 m6A regulators between patterns C1 and C2 was observed in the TCGA LIHC cohort (Fig. 5A and S5C). The expression of these m6A regulators (except for the eraser ALKBH3) in the pattern C2 showed higher expression level (Figure S5C). PCA also revealed that patients were well classified by the expression of 26 m6A regulators (Figure S5B). The Kaplan-Meier survival analysis showed that pattern C1 had a better prognosis than pattern C2 (Fig. 5B, log-rank test, P = 0.003). The ssGSEA was conducted in the TCGA LIHC cohort, and similar results, such as activated CD4+ T cells, activated CD8+ T cells, effector memory CD8+ T cells, and type 2 T helper cells, were observed (Figure S5D). The correlation between the immune characteristics of TME and the two m6A methylation patterns was then investigated further. We evaluated the association between two m6A methylation patterns and several immune infiltration-related signatures from the previous study [52] (Table S7). As shown in Fig. 5C, pattern C1 was enriched in CD8+ T effector, angiogenesis, cytolytic activity, and type II IFN response, etc. Pattern C2 had a higher expression of abundant epithelial-mesenchymal transition (EMT) 1 and EMT 3 than pattern C1. Furthermore, we constructed a heatmap with the expression of immune-related genes to visualize the relative expression level between patterns C1 and C2 (Fig. 5E). Ligands (such as CXCL10, CXCL9, VEGFA, and VEGFB), co-inhibitor CD276, and cell adhesions (ITGB2 and ICAM1) showed a relatively high expression in the pattern C2. In addiction, based on the two m6A methylation molecular subtypes of the TCGA LIHC cohort, we noticed that many patients with advanced HCC belonged to pattern C2, whereas the majority of patients with lower histology grades belonged to pattern C1, indicating that pattern C2 was involved in the progression of HCC (Fig. 5D and 5F). The stemness index based on mRNA expression (mRNAsi) is an indicator that describes how similar tumor cells are to stem cells on the basis of mRNA expression level [34]. Considering the strong connection of pattern C2 with later staging and histology grade subtype, we next investigated the correlation between m6A methylation molecular subtypes and mRNAsi. As expected, the mRNAsi was significantly evaluated in pattern C2, which further identified that pattern C2 may be involved in the malignant progression of HCC (Figure S6A).
The above findings suggested that HCC patients could be classified into two molecular subtypes based on the expression of 26 m6A regulators. In the two subtypes, pattern C1 represented a high proportion of CD8+ T cell infiltration and stromal cells with a favorable prognosis. Pattern C2 represented a high proportion of CD4+ T cells, tumor cell purity, and type 2 T helper cell infiltration with a poor prognosis.
Ma Phenotype-related Degs And The Construction Of The Ma Signature In Hepatocellular Carcinoma
The empirical Bayesian approach was used to identify DEGs between the two m6A modification patterns in order to further explore the molecular mechanisms underlying the two m6A methylation regulator patterns. A total of 2887 DEGs were identified as m6A-related DEGs (Figure S6B). GSEA between patterns C1 and C2 revealed that gene sets involved in cell cycle-related pathways, such as E2F targets, G2M checkpoint, MYC targets V1, TGF-β, and EMT, were significantly enriched in pattern C2, which was consistent with previous analyses (Fig. 6A). The Kaplan-Meier and univariate Cox regression analyses confirmed the prognostic value of 861 genes (Table S8). Unsupervised clustering analysis was performed on the expression of these 861 genes and classified HCC patients into two patterns, which we termed m6A methylation-related gene patterns (Fig. 6B and S6C). The clinical analysis revealed that patients with pattern B had a poor prognosis and a higher histology grade subtype (Fig. 6B). Kaplan-Meier survival analysis showed that pattern A had a better prognosis (Fig. 6C). Immune-related gene expression patterns in the two m6A DEGs-related gene patterns were remarkably similar to the two m6A modification patterns in the TCGA LIHC cohort (Fig. 6D). Angiogenesis and type I/II IFN response signatures were higher in gene pattern A patients, whereas EMT1/3, immune checkpoint, and MHC-II HLA signatures were higher in gene pattern B patients (Figure S6D). These results further supported the potential m6A modification patterns that there were two different molecular subtypes in HCC, which represented different clinical and immune characteristics. Subsequently, we established a m6A signature based on 861 m6A methylation prognosis-related genes. The regression coefficient of each gene was estimated applying LASSO Cox regression, and the optimal lambda (λ) was derived by 10-fold cross validation (Figure S7A and S7B). We divided the HCC patients into high- and low-risk groups using a median cutoff value of 0.8434 in the TCGA training set. In the IMvigor210 cohort, the optimal threshold value of 4.1580 was adopted to split patients into high-risk and low-risk groups. A Sankey diagram summarizes the relationship between the m6A methylation regulator cluster, the m6A gene pattern, risk subtype, and overall survival status (Fig. 7C; Table S9).
Prognostic Value Of The Ma Signature And Exploration Of Its Clinical Relevance
Survival study was performed to investigate the prognostic value of the m6A signature in HCC patients, and the results showed that the low-risk group had favorable overall survival in the TCGA cohort (the TCGA training set, P < 0.001; the TCGA test set, P < 0.001; the whole TCGA cohort, P < 0.001; Fig. 8A-8C), the ICGC cohort (P < 0.001; Fig. 8D), the gathered four HCC cohort (P < 0.001; Fig. 8E), and the GSE116174 cohort (P = 0.002; Fig. 8F). However, the GSE76472 cohort did not reach statistical significance (P = 0.541; Fig. 8H). In the IMvigor210 immunotherapy cohort, patients with a low-risk score exhibited markedly prolonged survival (P = 0.019, Fig. 8G). Growing evidence has demonstrated that tumor mutation burden (TMB) is a potential biomarker for predicting the response to checkpoint immunotherapy. Consequently, we investigated the correlation of the high- and low-risk groups with TMB in the TCGA LIHC cohort. As displayed in Fig. 8I, TMB was relatively higher in the high-risk group (P = 0.0016; Fig. 8I), indicating the relatively high level of neoantigen in the high-risk group. Moreover, mRNAsi was positively correlated with the risk score (R = 0.49, P < 0.001; Figure S7C), which might further explain why the high-risk group had a poorer prognosis. The risk score could also explain part of the clinical characteristics. The low-risk group had earlier staging and histology grade in the TCGA cohort (Fig. 7A, ***P < 0.001). The results of the ICGC cohort supported that of the TCGA cohort as well (Fig. 7B, ***P < 0.001). In both the TCGA and ICGC cohorts, a multivariate Cox regression analysis demonstrated that the m6A signature was an independent predictive factor for evaluating individual overall survival (Fig. 9A and 9B). The above results indicated that the m6A signature might be a prognostic biomarker that could effectively predict overall survival and part of the clinical features in HCC patients.
The Role Of The Ma Signatue In Predicting Immunotherapeutic Benefits
Potential biomarkers for predicting immunotherapeutic response mainly include PD-L1 expression, TMB, and microsatellite instability [21, 53]. Newly discovered predictors, such as TIDE and IPS, are increasingly recognized for evaluating immune responses. We also assessed the correlation of TIDE and IPS with risk score in the gathered HCC cohort and the TCGA LIHC cohort. The results indicated that the TIDE was substantially decreased in the low-risk group (the gathered HCC and TCGA cohort, both P < 0.001, Fig. 9C-9D, respectively) and IPS was significantly elevated in the low-risk group (the gathered HCC cohort and TCGA cohort, both P < 0.001, Fig. 9E-9F, respectively). These results suggesting that m6A methylation patterns might play a crucial role in regulating the immune response.