Analysis of m6A Regulators Expression Pattern in HCC
Expression pattern of 21 m6A-related genes (Table 1)was systematically analyzed in tumor specimen and paired normal sample from TCGA HCC cohort. We observed that the expression levels of most m6A regulators were significant distinct between tumor tissues and adjacent samples (Figure 1A and Figure 1B). Concretely, m6A-related genes including ALKBH5, EIF3A, FTO, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, KIAA1429, METTL16, METTL3, RBM15, RBM15B, WTAP, YTHDC1, YTHDC2, YTHDF1, YTHDF2, and YTHDF3 (all p <0.001) were dramatically higher in HCC specimen relative to adjacent normal liver samples. However, we didn’t find statistically significant distinction in terms of METTL14 as well as ZC3H13 (p = 0.06, 0.83, respectively). Furthermore, we analyzed the correlation between 21 m6A regulators to further elucidate the inherent relationship. Among these m6A regulators, the intrinsic connection between HNRNPC and HNRNPA2B1 is the most significant presented in Figure 1C. To further understand the m6A regulators expression patterns in the standpoint of biological procedures, we employed GO annotation on genes whose expression level abnormally regulated in HCC tissues. Figures 1D and 1E shown that upregulated m6A-related genes were mainly enriched in mRNA-related regulatory processes, including regulation of mRNA metabolic process and regulation of mRNA stability.
Mutation of m6A Regulators in HCC
The genetic alteration information of m6A regulators were explored employing TCGA HCC cohort on the cBioPortal database to uncover the potential influence of genetic alteration upon corresponding gene expression (Figure 1F). On the whole, we found that VIRMA had the highest alternation proportion and exhibited 9% genetic alteration, and the most common alteration manner was amplification.
Immune Cells Infiltration Subsets in Tumor Immune Environment of HCC
To assess the constitute of 22 TICs types, the CIBERSORT algorithm was employed in not only TCGA but also ICGC dataset. The overall fraction of immune cells in HCC was presented in Figure 2A and 2B. The highest proportion of TICs was resting CD4 memory T cells in TCGA cohort, whereas naive B cells accounted for the most abundant infiltrating immune cells in ICGC cohort, suggesting activated immune cells mediated in antitumor response may exert an opposing player in HCC tumorigenesis and progression. Figure 2C and 2D presented the distributions of 22 immune cell proportion together with HCC patients. To elucidate the clinical significance of TICs, we correlated components of 22 TICs with clinicopathological characteristics. We found that the distribution of resting dendritic cell had closely correlation with patient gender (Figure 2E). The constitute of regulatory T cells reduced significantly with advanced clinical stage (all P<0.05, Figures 2F), indicating regulatory T cells might serve as a suppressing role in HCC development. To estimate the prognostic predictive performance of TICs, we analyzed patient prognosis based on distinct TICs fraction. Taking advantage of Kaplan-Meier method, we found that activated NK cells was significantly correlated with better prognosis in TCGA cohort (Figure 2G, P=0.046). Likewise, activated NK cells had closely association with longer overall survival in ICGC cohort (Figure 2H, P=0.038). These results suggested that Tregs and activated NK cells may play nonnegligible roles in antitumor response of HCC.
Correlation of Consensus Clustering with Prognosis, Clinical Features and Tumor Immune Environment Characterization in HCC
To better reveal the clinicopathological significance of 21 m6A RNA methylation regulators, we clustered patients into two subtypes based on m6A regulators expression patterns. According to similarities displayed in m6A modulators, we observed that k = 2 had optimal clustering stability. An increasing trend of the cumulative distribution function (CDF) value was regarded as indicator of an excellent clustering (Figures S1A–C). To further support the result of the consensus clustering, principal-component analysis (PCA) was performed which presented cluster1/2 were non-overlapping and differentiated well (Figure 3A). Subsequently, we found overall survival of the cluster 1 were longer than cluster 2 using Kaplan–Meier analysis (Figure 3B, P=2.682e−04). Then, differences in clinicopathological variables between two subgroups were investigated. As a result, most m6A-related genes were remarkably downregulated in cluster 1 compared with cluster 2. Figure 3C presented that cluster 2 possessed significant correlation with female gender, and advanced clinicopathological stage (both p < 0.05). Therefore, the results of consensus clustering suggested the expression pattern of m6A modulators may act as key regulators in HCC malignancy.
To elucidate the correlation of m6A regulators with TIME of HCC, we analyzed the immune infiltration type and extent and calculate corresponding immunoscore of cluster 1/2. We explored whether there was distinction between two HCC subtypes regarding to the estimate score, immune score, stromal score and tumor purity. Our results presented that relative to cluster 1, cluster 2 obtained lower estimate, stromal and immune scores (Figure 3E) whereas higher tumor purity (Figure 3D). The difference of immune-related signature between cluster1/2 was that subtype 2 was closely correlate with higher aDCs and MHC class I, whereas subtype 1 had more B cells, Neutrophils, NK cells, pDCs, Cytolytic activity, DCs, Mast cells, TIL and Type I/II IFN Reponse, indicating that the distinction of m6A regulators expression pattern significantly correlated with TIME characterization of HCC (Figure 3F and Figure S1D). Figure S1E shown that each patient immune-related signature with corresponding immune scores in cluster1/2. Next, the differential subpopulation of infiltrating tumor immune cells between two different subtypes was identified. The results presented that cluster 1 had higher infiltration levels of B cells memory and Monocytes, while immune infiltration of follicular helper T cells and M0 Macrophages were remarkably lower (Figure 3G).
To investigate the involvement of m6A RNA methylation with ICB treatment, we
analyzed expression level of 47 ICB-related genes between two clusters. Compared with cluster 1 group, the majority of immune checkpoint blockade-related genes expression levels of cluster 2 were dramatically higher (i.e., PDCD1 and CTLA4, etc.; Figures 3H). Hence, the clustering results might contribute to reveal the complexity of TIME and predict ICB therapy outcome in HCC.
Construction of Prognostic Risk Signature
To further explore the prognostic predictive performance of m6A modulators in HCC patients, we conducted univariate Cox regression analysis on 21 m6A regulators expression levels based on TCGA database. We observed that 14 out of 21 m6A regulators had significant association with overall survival (Figure S2A). Among these 14 regulators, YTHDF2, YTHDF1, IGF2BP3, METTL3, RBM15B, HNRNPA2B1, KIAA1429, HNRNPC, WTAP, IGF2BP1, YTHDC1, RBM15 and IGF2BP2 were deemed unfavorable prognostic factors; whereas only ZC3H13 was regarded as a beneficial prognostic indicator. Then, LASSO algorithm was analyzed to screen the m6A RNA methylation regulators with the most powerful prognosis predictive ability (Figure S2B, S2C).
Finally, six m6A-related genes, including YTHDF1, YTHDF2, IGF2BP3, KIAA1429, METTL3, and ZC3H13, were identified to yield a m6A-based risk signature for HCC patients. Figure S2D presented the corresponding coefficients.
The risk score of HCC patients was calculated as following equation: risk score = (0.0262 ∗ expression level of YTHDF1) + (0.0577 ∗ expression level of YTHDF2) + (0.1192 ∗ expression level of IGF2BP3) + (0.027 ∗ expression level of KIAA1429) + (0.0795 ∗ expression level of METTL3) – (0.1018 ∗ expression level of ZC3H13).
Subsequently, each HCC patient obtained corresponding risk score and were randomized into high-/low-risk subgroups according to the median threshold.
Confirmation of Risk Prognostic Signature
Figure S3A displayed the distributions of six m6A regulators expression level with corresponding subgroups and patients. The allocations of risk score and dot pot of survival status in the TCGA HCC cohort suggested that high-risk HCC patients had shorter overall survival (Figures S3C, S3E). Besides, Kaplan–Meier curve corroborated that patients with low-risk possessed significant better prognosis than patients in high-risk group (P = 1.544e−04; Figure 4A). To assess the prognostic value of risk signature in HCC cohort, we performed ROC curve analysis. Area under curves of risk score signature at 3‐year survival times were 0.724, suggesting great sensitivity and specificity of the survival predictive ability (Figure 4C). Besides, results of univariate Cox model presented the hazard ratio (HR) of risk score was 3.713 (95% CI: 2.411−5.716; Figure S3G). And corresponding results were obtained in multivariate Cox regression analysis (HR = 3.386, 95% CI: 2.168−5.290; Figure S3I), suggesting risk score could act an independent indicator in HCC. Furthermore, the involvement of m6A-related genes with clinicopathological features was investigated and presented in the heatmap (Figure 4E). We observed that with advanced clinical stage(2 out of 6 P<0.05, Figures 4G) and advanced pathological grade (most P<0.05, Figures 4H), risk score significantly elevated.
Validation of Risk Prognostic Signature
To further estimate its external prognostic validity, we employed ICGC dataset (LIRI) as an external testing group. The ICGC-LIRI cohort with 231 HCC patients were classified into high-risk and low-risk subgroups according to the median threshold of TCGA dataset. The results presented the distributions of m6A regulators expression patterns, survival status, and risk score in the external validation cohort (Figure S3B, Figure S3D and S3F) and the combination set (Figure S4A, S4B and S4C). Same as the results from training set, Kaplan–Meier analysis presented that HCC patients with high-risk possessed significantly poorer prognosis relative to low-risk group patients in both validation cohort (Figure 4B, P = 6.69e−03) and the combination set (Figure S4D, P = 5.545e−05). The values of area under the ROC (AUC) was 0.76 in the external testing set (Figure 4D), suggesting good prognostic performance of risk prognostic signature among different populations. Consistent with results in the training group, risk signature as a prognostic factor independently affected overall survival in both univariable and multivariable regression analysis of not only validation group but also the whole cohort (Figure S3H, S3J, S4E and S4F). Subsequently, we plotted the heatmap to simultaneously presented clinical relevance (Figure 4F). The higher the risk score, the more serious the tumor stage (most P<0.05, Figures 4I).
Intimate Correlation of Prognostic Risk Score with Tumor Immune Environment Characterization of HCC
To further examine whether risk score can act as immune indicators in both TCGA and ICGC HCC datasets, we performed the correlation analysis of prognostic risk score with immune score (calculated using the ESTIMATE algorithm), ssGSEA signatures and TICs subtype and level (calculated via CIBERSORT method), and the expression level of 47 immune checkpoint blockade-related genes.
We found that low-risk patients obtained a higher stromal score compared with high-risk HCC patients in TCGA dataset but not ICGC cohort(Figure 4J, S5B). However, there was no significant difference regarding to immune score, estimate score and tumor purity (Figures S5A, S5B). Subsequently, we distinguished distinction of the immune-related signatures between these two subgroups. Combining two dataset ssGSEA results, the infiltration of aDCs, DCs, Th2 cells and some immune signatures like check-point, HLA molecule expression level, HLA molecule expression and MCH class I expression were significantly increased with elevated risk score (Figure 5A, 5B, S5C, S5D). Figure S5E and S5F shown that immune-related signature of each patient with corresponding immune scores in low-/high-risk groups in two datasets. The CIBERSORT algorithm results indicated that proportion of Tregs (regulatory T cells) was positively associated with risk score based on TCGA dataset(Figure 5C), whereas ICGC patients with high-risk presented lower gamma delta T cells, less M1 macrophages and more neutrophils relative to low-risk group (Figure 5D). Further correlation analysis presented that 25 of 47 (i.e., CTLA4, PDCD1, etc.,) immune check blockade-associated genes expression levels were significantly upregulated in patients with high-risk based on two datasets(Figure 5E, 5F). Above results indicated that m6A-based risk signature may provide a novel approach to elucidate the characteristics of immunity regulatory network and further forecast immunotherapy outcome in HCC.
Confirmation of Prognostic Value of m6A-based signature in HCC
Then, we analyzed a ROC curve and the value of AUC for the 1-, 2-, and 3-year OS were 0.746, 0.725, and 0.731, respectively, which indicating good predictive accuracy (Figure 6A). To explore whether m6A-based signature was the best prognostic indicator among various clinical characteristics, age, gender, stage and grade were extracted as the candidate prognosis predictive factors.
We integrated these clinical variables then conducted the AUC curve analysis for 1-, 2-, and 3-year OS and observed that risk signature obtained the most AUC (Figures 6B, 6C, and 6D). We generated a nomogram including risk score and stage to forecast prognosis of patients with HCC (Figure 6E). Age, gender and grade were rejected out of the nomogram because of their AUCs were less than 0.6. Calibrate curves indicated powerful prognostic predictive ability of 1-, 2- and 3-year OS in our nomogram plot (Figure 6F-H).
We performed a stratification analysis to validate whether m6A signature still had powerful prognostic predictive ability when HCC patients classified into various subgroups based on clinical characteristics. Relative to patients with low-risk, high-risk HCC patients presented poorer prognosis in both the early- and late-stage subgroups (Figures S6A, S6B). Similarly, m6A signature presented excellent prognostic prediction performance for patients in T1-2 or T3-4 status (Figures S6C and S6D), patients male gendered (Figures S6E), patients in grade 1-2 (Figures S6G), patients aged <=65 years (Figures S6I), patients in N0 status (Figures S6K), and patients in M0 status (Figures S6L). Meanwhile, we found the prognostic predictive ability of our signature was lost in patients in female gendered (Figures S6F), grade 3-4 (Figures S6H) or aged >65 years (Figures S6J). These results suggested that it can be an outstanding predictor in patients with HCC.
Functional Annotation of Prognostic Risk Signature
To further explore the potential player of m6A-based risk signature involved in HCC from a perspective of biological processes, we carried out GSEA in the low‐ as well as high‐risk subgroup. GSEA results shown that the high-risk score of signatures significantly enriched in pathways (i.e., non-small cell lung cancer, prostatic cancer, , the Wnt signal pathway, mTOR signal pathway, MAPK signal pathway, and the p53 signal pathway, etc., Figures S4G).
Correlation of Risk Signature with Immune Infiltration and ICB Key Molecules
Furthermore, we assessed whether this m6A-based signature correlated with infiltrating immune cells in TIME. The results shown that risk signature presented significant positive correlation with infiltrating B cells (r = 0.218; P = 2.752e−05), infiltrating CD4+T cells (r =0.200; P = 1.151e−04), infiltrating CD8+T cells (r = 0.209; P = 5.891e−05), infiltrating Dendritic cells (r = 0.305; P = 2.735e−09), infiltrating Macrophages (r = 0.404; P = 8.609e−16) and infiltrating Neutrophils (r = 0.349; P = 6.339e−12; Figures 6A-D). These results contribute strong evidence to validate that m6A-based risk signature had close relationship with infiltrating immune cell type and level in HCC.
Subsequently, we correlated six key immune checkpoint inhibitors genes(PDCD1, CD274, PDCD1LG2, CTLA‐4, HAVCR2, and IDO1) (38-40). And we analyzed the correlation between ICB key targets and m6A-based signature to reveal the potential player of risk signature in the ICB treatment of HCC (Figure 6G). We found that m6A risk signature was significantly positive correlated to CTLA4 (r= 0.15; P= 0.0013), HAVCR2 (r= 0.19; P= 5.2e−05), IDO1 (r= 0.093; P= 0.05), PDCD1 (r= 0.11; P= 0.021) and PDCD1LG2 (r = 0.12; P = 0.0097; Figures 6H-L), suggesting m6A signature might act an nonnegligible role in the prediction of responsiveness to ICB treatment in patients with HCC.
ZC3H13 Independently Affected Prognosis and Correlates with Immune Cell Infiltration ICB key Genes
ZC3H13 was only one m6A regulator whose expression level was downregulated among the prognostic m6A-related genes. Therefore, the role of ZC3H13 in HCC was explored in further experimental validation. We compared ZC3H13 expression level between normal tissues and tumor samples based on TCGA and GTEx data. Relative to tumor tissues, ZC3H13 expression level was lower in adjacent normal specimens (Figure 7A). Taking advantage of qRT-PCR, we determined ZC3H13 expression level in four distinct HCC cell lines and human hepatic cell line. Consistent of results of online database, ZC3H13 was downregulated in cancer cells relative to normal cell (Figure 7B). To further assess the prognostic value of ZC3H13 in HCC, Kaplan–Meier analysis were conducted between ZC3H13 low- and high-expressed patients. As presented in Figure 7C, higher ZC3H13 expression level significantly suggested longer overall survival time (P = 2.514e−06).
To elucidate the intrinsic relationships between infiltrating immune cell and ZC3H13 expression level, we analyzed correlation between ZC3H13 expression and immune cell infiltration level via using TIMER. We found that expression of ZC3H13 presented significant correlation with CD4+ T cells (r= 0.125; P = 1.97e−02), CD8+ T cells (r= 0.171; P = 1.47e−03), Myeloid Dendritic cells (r= 0.124; P = 2.11e−02), Macrophages (r= 0.134; P = 1.30e−02), and Neutrophils (r= 0.244; P = 4.62e−06; Figures 7D-H).
Then we performed the correlation between the ZC3H13 and ICB key targets adjusted by tumor purity by TIMER to investigate the potential player of ZC3H13 in ICB treatment of HCC. TIMER results presented ZC3H13 was significantly positive correlated to CD274 (r= 0.437; P = 1.73e−17), HAVCR2 (r= 0.14; P = 9.34e−03), IDO1(r= 0.113; P = 3.57e−02) and PDCD1LG2 (r = 0.187; P = 4.90e−04; Figures 7I-L), suggesting ZC3H13 may exert a vital player in ICB treatment of HCC.
Close Connection between ZC3H13 and Tumor Immune Microenvironment Characterization
To further elucidate the relationship between ZC3H13 and TIME characteristics in HCC, we analyzed the correlation of ZC3H13 expression value with immune scores and tumor purity (emplying the ESTIMATE method), ssGSEA signatures (using GSEABase algorithm), TICs subpopulation and level (via CIBERSORT tool) and the expression level of 47 immune checkpoint blockade-related genes. HCC patients were classified into high-/low- ZC3H13 subtypes based on the median ZC3H13 expression level. ESTIMATE results indicated that patients with lower ZC3H13 expression had a significant lower stromal score relative to patients in high-ZC3H13 group in TCGA cohort but not ICGC dataset. In terms of estimate score, immune score and tumor purity, however, there was no remarkable distinction between these two groups (Figures 9A, 9B). Subsequently, we identified difference of enrichment in immune-related signatures between two different subgroups. Based on two dataset ssGSEA outcomes, the infiltration fraction of Th2 cells, T cell co-inhibition, and check-point were significantly increased when risk score declining, whereas IFN-response type-II was positively correlated with ZC3H13 expression level (Figure 9C, 9D). The CIBERSORT analysis results of TCGA cohort shown that the proportion of activated NK cells was significantly higher in patients with low-risk (Figure 5E). However, there was no significant difference in ICGC dataset (Figure 5E). Taking advantage of correlation analysis, we found that 3 immune check blockade-related genes (i.e., TNFSF14, TNFRSF4, KIR3DL1) were significantly upregulated, but TNFRSF14 and TNFRSF18 were lower, in the high-risk group based on two datasets(Figure 9G, 9H). These results pointed out that ZC3H13 may serve as a key indicator in TIME characterization and immunological reaction in HCC.