Identification of m6A regulators in hepatocellular carcinoma.
The 13 m6A regulators expression in the TCGA database is shown in Figure 1A and 1B. Except for METTL14 and ZC3H13, the expression of the other 11 genes was strikingly different among tumor group and normal group. The expression of ALKBH5 (p=0.001), WTAP (p<0.05), YTHDF1, YTHDF2, FTO, HNRNPC, YTHDC1, YTHDC2, RBM15, KIAA1429, and METTL3 (P<0.001) in HCC tissues was upregulated compared with normal liver tissues (Figure 1A, 1B). 13 m6A regulators expression were positively correlated with each other (Figure 1C).
Consistent clustering of m6A regulators determined two clusters
In this study, we used the Consensus Cluster Plus package to extract 50 normal liver tissues and 374 liver cancer tissues. On the basis of the expression similarity of m6A regulators in the TCGA dataset, k = 2 seemed to have the smaller Cumulative distribution function (CDF) value (Figure 2 B, C). We divided the patients into two groups (Figure 2A). Next, we used PCA to analyze the two clusters. The results showed that cluster 1 and cluster 2 were distinct. (Figure 2D)
Correlation between consensus cluster classification and clinical characteristics
We analyzed the OS curves of two groups of patients. Kaplan-Meier survival curve was used to analyze OS rate of HCC patients between two clusters. In this study, we found that the OS of cluster 1 was striking shorter than cluster 2 (Figure 3A). Most m6A regulators were highly expressed in the cluster 2. Patients of HCC in the cluster 1 had higher scores compared with the cluster 2. In this study, we investigated the positive relationship of m6A regulator and the pathological characteristics such as age, gender, grade, stage, TMN status of HCC. (Figure 3B)
Risk Signature and Poor prognosis of m6A Regulators
Univariate Cox regression analysis was used to the expression of the TCGA dataset. The results showed that HNRNPC (HR= 1.024, 95% CI= 1.009–1.039), RBM15 (HR= 1.264, 95% CI = 1.031–1.548), WTAP (HR = 1.086, 95% CI = 1.029–1.147), KIAA1429 (HR = 1.158, 95% CI = 1.070–1.253), YTHDF2 (HR = 1.117, 95% CI = 1.070–1.165), METTL3 (HR = 1.252, 95% CI = 1.128–1.390), YTHDC1 (HR = 1.091, 95% CI = 1.015–1.173) and YTHDF1 (HR = 1.080, 95% CI = 1.045–1.115) corresponded to lower survival rate in patients. Instead, increased expression of ZC3H13 (HR = 0.907, 95% CI = 0.833–0.988) was associated with a higher survival rate in patients (Figure 4A). The 13 genes in the TCGA dataset were analyzed using the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm. (Figure 4B, C). We separated the patients into a low-risk group and a high-risk group. The results indicated that the high-risk group of patients had poor survival (Figure 4D).
Correlation between prognostic risk score and clinicopathological characteristics
We investigated the role of the five genes in the high and low-risk patients, and the pathological characteristics including age, stage, TMN status of liver cancer (Figure 5A). Higher expression of KIAA1429, YTHDF1, YTHDF2, and METTL3 and lower expression of ZC3H13 in hepatocellular carcinoma (Figure 5A). Using the ROC curve, we predicted prognostic risk scores and 3-year survival rates in HCC patients. The consequences showed that risk score predicted 3-year survival in patients (AUC = 0.619) (Figure 5B). Univariate and multivariate Cox regression analysis confirmed that the risk score, stage, T status were all related to OS (Figure 5C,D). The clinicopathological characteristics of primary liver cancer, and KIAA1429, YTHDF1, YTHDF2, and METTL3 were risk factors for poor prognosis of liver cancer patients.
Relationship between KIAA1429, YTHDF1, YTHDF2, and METTL3 expression and tumor-infiltrating immune cells
Independent tumor-infiltrating lymphocytes play a pivotal role in predicting prognosis (32). This study used TIMER to analyze the correlation between KIAA1429, YTHDF1, YTHDF2, and METTL3 expression and the level of immune invasion of liver cancer. KIAA1429 expression was associated with B cells (p-value = 1.39 × 10−7), CD4 + T cells (p-value =4.20 ×10−7), Neutrophil (p-value = 2.46 × 10 −7), and Dendritic cells (p-value =1.85 × 10−7). YTHDF1 expression was associated with B cells (p-value = 2.70×10−11), CD4 + T cells (p-value =5.11×10−16), Macrophage (p-value = 5.23×10 −21), Neutrophil(p-value = 2.25×10 −18), and Dendritic cells(p-value =1.36 × 10−14). YTHDF2 expression was associated with CD4 + T cells (p-value =8.97 × 10−13), Macrophage (p-value = 7.23 × 10 −16), Neutrophil (p-value = 3.02 ×10 −21), and Dendritic cells (p-value =6.71 × 10 −11). METTL3 expression was associated with B cells (p-value = 9.51×10−7), CD4 + T cells (p-value =5.73 ×10−14), Macrophage (p-value = 1.25×10 −11), Neutrophil(p-value = 8.70×10 −9),and Dendritic cells(p-value =6.21×10−9) (Figure 6A-6D). We further investigated whether there is a difference in the tumor immune microenvironment between high and low levels of the four m6A regulators in HCC. The CIBERSORT algorithm was applied to 22 immune cell subtypes to evaluate the difference in their expression levels between the high expression group and the low expression group (Figure 7A-7D). There was no significant difference in immune infiltrating cells between the KIAA1429 high expression group and the low expression group. Compared with the low expression group, the neutrophils in the YTHDF1 high expression group increased (p-value=0.039), the neutrophils in the YTHDF2 high expression group increased (p-value = 0.023), and the METTL3 memory B cells increased (p-value=0.006). The results indicate that KIAA1429, YTHDF1, YTHDF2, and METTL3 may play a crucial part in immune invasion of liver cancer.