In this study, we used LASSO (Least Absolute Shrinkage and Selection Operator) regression to screen out 7 genes (ALKBH5, FTO, HNRNPA2B1, IGF2BP2, IGF2BP3, YTHDF1, YTHDF3) significantly associated with prognosis and constructed a risk scoring model based on regression coefficients (Figure 2A, Figure 2B). Then, the risk scoring model was validated using cox proportional-hazards model (Figure 2C). The forest plot showed the effect of the references, which were age, gender and risk score, on the outcome which defined as death. We concluded that gender, in the TCGA-LAML cohort, has no effect on the survival of patients (p>0.05). Both age and risk score were factors that affect survival (p<0.001). Furthermore, our risk scoring model has an ideal area under the curve (AUC), which is 0.752 in 1-year survival, 0.789 in 3-year survival, and 0.864 in 5-year survival (Figure 2D). The Kaplan-Meier curve based on the risk scoring model also showed statistical significance (p<0.0001). Among TCGA-LAML patients, the prognosis of the high-risk group was significantly worse than that of the low-risk group (Figure 2E), indicating the relationship between the regulators and the prognosis. In addition, the GSE37642 cohort validated the risk scoring model (Figure 2F), and in GSE37642, the high-risk and low-risk groups were also significantly separated, showing statistical significance (p<0.01).
Differential expression analysis of m6A regulators in AML
We obtained the expression data of m6A from the TCGA and GTEX databases, using TCGA-LAML as the tumor group and GTEx as the normal control group. Figure 3 shows that the differential expression level of m6A related regulators between tumor group and normal group. As we can see, the expression levels of METTL3, MTEEL14, WTAP, KIAA1429, RBM15, RBM15B, ZC3H13, METTL16, FTO, ALKBH5, YTHDF3, YTHDF2, YTHDF1, YTHDC1, HNPRNPC, HNRNPA2B1, IGF2BP1, IGF2BP2, and IGF2BP3 were significantly different in AML and normal tissues (p < 0.05). Figure 4 showed that the differential expression heatmap of 19 regulators of m6A which indicated that the HNRNPA2B1, WTAP, HNPRNPC, IGF2BP3 and KIAA1429 were highly expressed in normal tissues compared with the tumor group. It also showed that the expression levels of IGF2BP2, RBM15, RBM15B, YTHDF3, YTHDF2, YTHDF1, YTHDC1, ZC3H13, MTEEL14, METTL16, METTL3, FTO, ALKBH5 were higher in AML samples than normal group. The overexpression or reduction of IGF2BP1 was unclear between the two groups.
Function Enrichment Analysis and Protein-Protein Interaction Network
The Gene Ontology molecular function enrichment (GO:MF) analysis indicated that the m6A-associated regulators were enriched in RNA binding, translation regulator activity, methyltransferase activity, transferase activity, RNA splicing, nucleotide binding, S_adenosylmethionine_dependent methyltransferase activity, oxidoreductase activity, and telomeric DNA binding, which showed that m6A regulators were involved in the metabolic process of RNA (Figure 5A). The enrichment of these functions may be the underlying mechanism of the pathogenesis of AML. Figure 5B showed the protein-protein interaction network of m6A regulators. As we can see, the proteins of m6A regulators interacted with each other, not only are they independent. It should be noted that they formed part of the cellular response network and regulated cell behavior to a certain extent.
Figure 6 presented the CNV analysis of m6A related regulators in AML. The pie plot summarized the CNV of m6A regulators’ genes in the AML, in which pointed out HNRNPA2B1 had 8.9005236% of copy number deletion in total, including heterozygous and homozygous deletion (Supplement Table1). Similarly, IGF2BP3, ALKBH3, RBM15B, FMR1, RBMX, FTO and IGF2BP1 also showed copy number deletion. On the other hand, copy amplification also includes heterozygous amplification and homozygous amplification. Through Figure 6A, we can see that the heterozygous amplification of YTHDF3 was very significant accounted for a percentage of 12.565445 and ZC3H13, YTHDF2, YTHDF1, RBM15, METTL14, HNRNPC, YTHDC1, WTAP, METTL3, IGF2 IGF2BP2, FTO, RBMX, FMR1 also showed copy amplification. It is worth noting that in AML, METL16 and ALKBH5 only had copy deletion, but no copy amplification and RBM15, WTAP, YTHDF1, YTHDF3 only had copy amplification, but no copy deletion (Figure 6B and Figure 6C). Figure 6D showed the relationship between CNV and mRNA expression through correlation analysis which was not very obvious, because the Spearman correlation coefficients were all less than 0.03 (Supplement Table). Furthermore, Figure 6E provided that the LogRank tests were performed to test the difference in AML. We can see that the CNV of METTL3, IGF2BP2, ZC3H13, RBM15B, METTL14, FTO, IGF2BP1, YTHDC1, YTHDF2 genes and the overall survival of AML were statistically different (p <0.05).
The results of correlation analysis between methylation level and mRNA expression level were shown in Figure 7A. We found that the methylation of IGF2BP2, FTO, ALKBH5, IGF2BP3 is negatively correlated with the expression level, and the correlation between the methylation of other regulators and the expression level is not very obvious. According to Figure 7B, there is a statistical significance between the methylation of ZC3H13, IGF2BP1, and HNRNPC and survival. Their hazard ratios are all less than 1, which means lower methylation has higher risk of death.
Pathway Activity Analysis
We conducted the difference of genes expression between pathway activity groups (activation and inhibition) that defined by pathway scores through GSCALite. Global percentage (Figure 8A) of cancers in which a gene has effect on the pathway among 32 cancers types, presents the percentage of cancer type that are activated or inhibited. We found that m6A related regulators can activate or inhibit apoptosis, cell cycle, DNA damage response, epithelial-mesenchymal transition (EMT), hormone AR, hormone ER, PI3K/AKT, RAS/MAPK, RTK, TSC/mTOR pathways. Figure 8B shows the m6A regulators that have function in at least 5 cancer types. Through the pan-cancer analysis, we can see that in the apoptosis pathway, FTO, METTL14 and YTHDC1 mainly play an inhibitory role and HNRNPA2B1, HNRNPC, IGF2BP3, RBM15 and WTAP mainly play a role in activation. In the cell cycle pathway, FTO strongly inhibit it, yet HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP3, KIAA1429, METTL3, RBM15, RBMX and YTHDF1 mostly activate it. In DNA damage response, HNRNPA2B1, HNRNPC, IGF2BP3, RBM15B and RBMX mainly play a role in activation. In epithelial-mesenchymal transition, HNRNPC plays a role in inhibition, and FTO and IGF2BP2 mainly play a role in activation. In hormone AR or ER pathways, HNRNPA2B1 and KIAA1429 activate the hormone AR, yet inhibit the hormone ER. Other m6A regulators are compatible with the regulation of estrogen and androgen receptors. In PI3K/AKT pathway, METTL14 activates it while IGF2BP2 inhibits it. In RAS/MAPK pathway, HNRNPA2B1, HNRNPC and RBMX play an inhibitory role and none of the m6A related regulators evidently activate the process of RAS/MAPK. In RTK pathway, HNRNPC plays an inhibitory role while FTO, METTL14 and YTHDF3 activate it. And we can see that YTHDF2 mainly activate the pathway of TSC/mTOR. In all, it should be noted that the involvement of m6A in the tumor process remains multi-pathway and multi mechanism.
Drug Sensitivity Analysis
As we can see, Figure 9A summarizes the correlation between gene expression and the sensitivity of Genomic of Drug Sensitivity in Cancer (GDSC) top 16 drugs. The expression of IGF2BP2 had a negative relationship with the sensitivity of trametinib and docetaxel and a positive relationship with the sensitivity of Navitoclax, KIN001-102, NPK76-II-72-1, I-BET-762 (FDR <0.05) (Supplement Table2). Moreover, Figure 9B and Supplement Table3 presents the correlation between m6A related gene expression and the sensitivity of Genomics of Therapeutics Response Portal (CTRP) top 30 drugs. The higher the expression of IGF2BP2, the strong the drug sensitivity of LRRK2-IN-1, serdemetan, QW-BI-011, PRIMA-1, BRD-A94377914, belinostat, BIX-01294, PX-12, Panobinostat (FDR <0.05). It is important to point out that m6A related genes may become a new target for tumor therapy.
In this study, we analyzed the co-expression within m6A regulators and the co-expression between m6A and 32 common mutations in AML (Figure 10). First, in terms of the co-expression of m6A regulators, we found that METL14 has a moderately positive correlation with KIAA1429, ZC3H13, and YTHDC1. KIAA1429 has a moderately positive correlation with ZC3H13 and YTHDC1. METTL3 is negatively correlated with the expression of ALKBH5, and ALKBH5 is negatively correlated with the expression of YTHDC1. Then, in the co-expression between the common mutations of AML and m6A regulators, we can see that METTL3 is positively correlated with ASXL1 and SF3B1; METL14 has a positive correlation with NF1, PPM1D; KIAA1429 shows positive correlation with NF1, NRAS and STAG2, and negative correlation with CEBPA; there exists a correlation between RBM15B and IDH2; ZC3H13 has a correlation with NF1 and STAG2; METL16 has a correlation with NF1 and NPM1. There is a correlation between ALKBH5 and TP53, and ALKBH5 remains negatively correlated with SF3B1; YTHDF1 is negatively correlated with ZRSR2; YTHDF2 is correlated with NPM1; YTHDF3 is negatively correlated with IDH2, and is positively correlated with NRAS; YTHDC1 is positively correlated with NF1, STAG2; HNRNPC is positively correlated with CSF3R and negatively correlated and positively correlated with NPM1; HNRNPA2B1 is negatively correlated with BCORL1 and MPL; IGF2BP3 is negatively correlated with CALR and CEBPA and positively correlated with SETBP1.