Global analysis of N6-methyladenosine functions and its disease association using deep learning and network-based methods

Abstract

N6-methyladenosine (m6A) is the most abundant methylation, existing in >25% of human mRNAs. Exciting recent discoveries indicate the close involvement of m6A in regulating many different aspects of mRNA metabolism and diseases like cancer. However, our current knowledge about how m6A levels are controlled and whether and how regulation of m6A levels of a specific gene can play a role in cancer and other diseases is mostly elusive. We propose in this paper a computational scheme for predicting m6A-regulated genes and m6A-associated disease, which includes Deep-m6A, the first model for detecting condition-specific m6A sites from MeRIP-Seq data with a single base resolution using deep learning and Hot-m6A, a new network-based pipeline that prioritizes functional significant m6A genes and its associated diseases using the Protein-Protein Interaction (PPI) and gene-disease heterogeneous networks. We applied Deep-m6A and this pipeline to 75 MeRIP-seq human samples, which produced a compact set of 709 functionally significant m6A-regulated genes and nine functionally enriched subnetworks. The functional enrichment analysis of these genes and networks reveal that m6A targets key genes of many critical biological processes including transcription, cell organization and transport, and cell proliferation and cancer-related pathways such as Wnt pathway. The m6A-associated disease analysis prioritized five significantly associated diseases including leukemia and renal cell carcinoma. These results demonstrate the power of our proposed computational scheme and provide new leads for understanding m6A regulatory functions and its roles in diseases.

Fig 1. Flowchart of our proposed prediction pipeline.

Fig 2. Performance of Deep-m⁶A.

(A) shows the IP reads count coverage in the 101-nt positive and negative training samples. The x-axis denotes the relative location of the nucleotides in the sample sequences, where the 51st position is the center A of the DRACH motif and the 1^st position is the 5’ end of the sequence. The IP expression level represents the IP reads input RC_norm. (B) and (C) are the ROC and PR curves of Deep-m⁶A-S and Deep-m⁶A obtained from the 10-fold CV on the HEK293 training data. (D) and (E) are performances of Deep-m⁶A, Deep-m⁶A-S, SRAMP-Mature and SRAMP-Full model on the independent MOLM13 data. The number after the method names are corresponding area under the curve (AUC).

Fig 3. Analysis of predicted m⁶A sites in 75 human samples.

(A) Distribution of the number of occurrence of a predicted single-base m⁶A site in the 75 samples. (B) Distributions of all predicted m⁶A sites and sites appeared in more than 12 samples in a meta-mRNA. (C) Distributions of the absolute revised Fisher’s z transformed Pearson correlations of the 49 consensus m⁶A-regulated genes, 709 m⁶A-regulated genes identified by HotNet2, and other remaining candidate genes. (D) GO BP and KEGG pathway enrichment result for all 709 m⁶A-regulated genes. Gene count means the number of genes involved in the corresponding terms and P_Benjamini is the adjusted FDR of the enrichment p-value.

Fig 4. The architecture of the proposed CNN model for Deep-m⁶A.

The input matrix is M_sr for Deep-m⁶A and M_s for Deep-m⁶A-S.

Fig 5. Gene-disease heterogeneous network.

The top network is gene-gene interaction network, the bottom network is disease-disease similarity network and they are connected by gene-disease relationship (dashed grey lines). The orange gene nodes denoted the m⁶A regulated genes and the green disease nodes are m⁶A regulated genes correlated diseases.