As the most prevalent mammalian mRNA epigenetic modification, N6-methyladenosine (m6A) has been shown to possess important post-transcriptional regulatory functions. However, the regulatory mechanisms and functional circuits of m6A are still largely elusive. To help unveil the regulatory circuitry mediated by mRNA m6A methylation, we develop here m6A-Driver, an algorithm for predicting m6A-driven genes and associated networks, whose functional interactions are likely to be actively modulated by m6A methylation under a specific condition. Specifically, m6A-Driver integrates the PPI network and the predicted differential m6A methylation sites from methylated RNA immunoprecipitation sequencing (MeRIP-Seq) data using a Random Walk with Restart (RWR) algorithm and then builds a consensus m6A-driven network of m6A-driven genes. To evaluate the performance, we applied m6A-Driver to build the context-specific m6A-driven networks for 4 known m6A (de) methylases, i. e. , FTO, METTL3, METTL14 and WTAP. Our results suggest that m6A-Driver can robustly and efficiently identify m6A-driven genes that are functionally more enriched and associated with higher degree of differential expression than differential m6A methylated genes. Pathway analysis of the constructed context-specific m6A-driven gene networks further revealed the regulatory circuitry underlying the dynamic interplays between the methyltransferases and demethylase at the epitranscriptomic layer of gene regulation.
Methylation, as a significant epigenetic modification of nucleic acids, regulates gene expression, influences grows and development of plants and animals, and is closely related to the occurrence and development of disease. The epigenetic regulatory mechanisms and physiological functions of DNA methylation have been well established through intensive studies in simple model organisms to human in the past decade [1–3]. However, RNA methylation, even though prevalent in many organisms, has long been considered to have little functional relevance. The discovery of obesity-associated FTO as a demethylase [4] of mRNA N6-methyladenosine (m6A) revealed that mRNA m6A methylation can be reversed and is thus a highly dynamic phenomenon. This discovery sparked the surged interests in study the prevalence of m6A in different cells and the functions of m6A. Subsequently, using methylated RNA immunoprecipitation sequencing (MeRIP-seq) technique [5–7], transcriptome-wide distribution of m6A in mammalian cells was profiled [6,7], revealing for the first time a widespread occurrence of m6A in >25% transcripts. m6A was also shown to be enriched around the stop codon of RNA transcripts and conserved between people and mouse [6,7], implicating a potential role played by m6A in post-transcriptional regulation [6,8, 9]. Since then, m6A has been shown to have a number of important biological functions, including promoting RNA degradation [10], regulating RNA stability by modulating binding of RNA binding proteins [6,11,12], and controlling translation efficiency [13–17]. Meanwhile, the identification of m6A methyltransferases and demethylases [4,18–20] further revealed the regulators of epitranscriptome. We now know that the m6A methyltransferase complex consists of METTL3, METTL14, and WTAP and functions as m6A" writers" in eukaryotes [9,18,21]. In contrast, FTO and ALKBH5 are identified to be de-methyltransferase, or m6A" erasers" [4,9, 22], indicating that mRNA m6A methylation is a dynamic process [4] and directly regulated by a number of methylases and demethylases [23]. Knockdown studies of these (de) methylases further revealed their involvement in many significant physiological processes including obesity [24–26], synaptic signaling [27], cancer [28,29], sperm development [22], stem cell differentiation [30], circadian periods [31], yeast meiosis [32,33], and stem cell pluripotency [34–36]. Although these studies together greatly improve our understanding of the reversible mRNA m6A methylation, the regulatory mechanisms and functional circuitry of m6A are still largely elusive. Currently, MeRIP-Seq is the most widely adopted high throughput approach for measuring transcriptome-wide m6A methylation [6,7, 37]. To obtain a transcriptome-wide m6A profile, MeRIP-Seq produces two sets of samples, i. e. , IP and input samples. While IP samples include sequencing reads from m6A methylated RNA fragments pulled down with anti-m6A antibody, input samples measure the basal abundance of all transcripts, which are used as background for assessing the enrichment of methylated fragment. Detecting m6A methylation site or" peak detection" from MeRIP-Seq data can be achieved by comparing the enrichment of reads in the IP samples vs. those in the input samples. Several algorithms including exomePeak have been developed for m6A peak detection [38–41]. After the methylation sites are identified, differential m6A methylation (DmM) analysis can be also performed in a case-control study to further identify the dynamic m6A sites whose methylation status is significantly different under two experimental conditions. Algorithms such as exomePeak [42] and MeTDiff [43] have also been developed for this purpose. While peak detection and DmM analysis are essential steps for m6A bioinformatics analysis, they do not yet provide direct information about the functional relevance of m6A. We focus in this paper on predicting m6A-driven genes (mDrGenes) and the m6A-driven gene interaction network (mDrNet). Specifically, we refer mDrGenes as genes whose mRNAs harbor DmM sites or differential m6A methylation genes (DmMGs), and thus may be under dynamic epitranscriptomic regulation and be functionally significant to the biological context of interest. Conceivably, when data is available, mDrGenes can be conveniently identified by first predicting the DmMGs and then assessing their functional significance by using functional networks such as Protein-Protein Interactions (PPI) network or biological pathways. However, identifying functionally significant DmMGs when there are replicates can be nontrivial. The challenge arises as a result of technical and biological bias, where significant DmMGs identified in some replicates might not be significant in other replicates. Existing algorithms for DmM analysis such as exomePeak and MeTDiff all devise different methods ranging from taking consensus DmM sites [42] to statistically modeling of replicate samples [43] to mitigate this bias. While they can help detect robust DmM sites, these DmM sites might not be functionally significant DmMGs. As our goal emphasizes on detecting functional significance, an approach that can address this bias in assessing functional significance is more desirable and likely to better identify the m6A-driven genes and network. To address the aforementioned issue, we propose in this paper m6A-Driver, an algorithm that predicts mDrGenes by evaluating the consistency of RNA differential methylation from a functional network perspective. Specifically, rather than predicting DmMGs directly, m6A-Driver first performs DmM analysis on every possible replicate set (RS) independently, where each RS includes two IP-input pairs, one from the treated/case condition and the other from the untreated/control condition. Then, a DmM functional network is constructed for each RS by searching the significant interactions with DmMGs in PPI network using a Random Walk with Restart (RWR) algorithm. We adopt PPI network here to model functional interactions of m6A mediating genes because m6A has been shown to regulate the process of translation [13–17], in addition to its influence on gene expression. Finally, a consensus m6A-driven gene network is built by taking all the significant reoccurring interactions. By assessing the consensus among RS networks as opposed to RS DmMGs, m6A-Driver effectively addresses the sample bias that impacts functional prediction. m6A-Driver was applied to four case-control studies that investigate the functions of the component of methyltransferase complex (METTL3, METTL14, and WTAP) and demethylase (FTO). In the end, m6A-driven gene networks were constructed for each (de) methylase together with an integrated network for the complete m6A methyltransferase complex. We showed that the predicted m6A-driven genes have higher degree of differential expression and more explicit functional relevance than DmMGs identified directly by previous approaches. These results demonstrate the effectiveness of m6A-Driver in prioritizing functional significant m6A-driven genes from m6A sequencing data.
The algorithm of m6A-Driver consists of four steps, depicted in Fig 1, with the first three steps implemented in each RS and step 4 performed to combine the results from all RSs. In step 1, exomePeak [42] is applied to detect DmMGs in each RS. In step 2, for each RS, the Random Walk with Restart (RWR) is performed using every DmMG as the seed node separately to search for their closely interacting genes in the PPI network. In step 3, the topological and biological significance of these DmM interacting genes are assessed and the genes that are determined to be insignificant are filtered out. The topological significance is estimated by their occurrence as top nodes prioritized by the same RWR algorithm in 100 random networks generated with the same topological structures. Meanwhile, the biological significance is evaluated by the length of their shortest path to the initial node (the seed DmMG). RS-specific DmM interacting networks, each consisting of significant interacting genes, are constructed at the end of step 3. Finally, in step 4, an mDrNet is constructed by assessing the interaction recurrence across all RSs. The genes that make up the nodes of the mDrNet are predicted mDrGenes determined by exomePeak. In this way, we extract a set of mDrGenes, or functionally relevant genes driven by m6A and a network that depicts the functional relationship of mDrGenes. We applied m6A-Driver on 4 MeRIP-seq datasets, i. e. , FTO knockdown dataset (KD-FTO) [23], METTL3 knockdown dataset (KD-METTL3), METTL14 knockdown dataset (KD-METTL4), and WTAP knockdown dataset (KD-WTAP) [44]. KD-FTO dataset is obtained from [27] that profiles m6A in FTO gene knockdown mice and their wild-type littermate. There are 12 samples (3 IP replicates paired with 3 input replicates for FTO knockdown mice and 3 IP replicates paired with 3 input replicates under wild-type littermate). It was divided into 9 sets of biological replicates and each biological replicate set (RS) contains two IP samples respectively from a FTO knockdown mouse and a wild-type (WT) littermate and two corresponding input samples from the two mice. KD-METTL3, KD-METTL14 and KD-WTAP datasets are from a recent study, which shows that m6A regulates mRNA stability [44]. Each dataset contains 8 samples, 2 IP replicates paired with 2 input replicates from the knockdown HeLa cells and 2 IP replicates paired with 2 input replicates from untreated HeLa cells. Similar to KD-FTO, samples in each of the three datasets are then divided into 4 RSs, each of which contains two IP samples from the knockdown HeLa cells and untreated HeLa cells respectively and two corresponding input samples. We first predicted the DmM sites in each dataset using exomePeak. As the technical limitation of MeRIP-Seq can lead to high sample bias, making the prediction results less reliable, we then set out to check the quality of the prediction results. First of all, the specificity of the predictions by exomePeak and MeTDiff on these datasets has been evaluated in a previously published paper [43], which shows that the false positive rates for all these datasets can be controlled and there are high specific DmM sites predicted in all these datasets. Next, we further examined the predictions of the three m6A methylase knockdown datasets, where we created a set of pseudo control and pseudo knockdown sequencing samples by scrambling the samples of a dataset so that the pseudo control and knockdown samples are both made up by a real control replicate and a real knockdown replicate. We then performed exomePeak on both the real dataset and the pseudo dataset and examined the prediction specificity at different thresholds by comparing the ratio of predicted DmM sites (or reported true positive rate, RTPR) in the real dataset and those in the pseudo dataset (or false positive rate, FPR) using a ROC-like curve. As is shown in S1 Fig, the percentage of DmM sites in the real datasets is much higher than that in the pseudo dataset at different thresholds in all the 3 datasets. Taken together, these results demonstrate that the false positive rates in these datasets can be controlled and the exomePeak prediction results are of good specificity for subsequent analysis. The reference network, PPI network, is built from the most recent version of PPI data from BioGRID (release 3. 4. 128, compiled on August 25th, 2015) [45]. Based on the binary interactions, we removed the isolated proteins and self-interaction proteins to establish a PPI network with a total of 16,062 proteins and 152,676 interactions. Jia et al. have proposed the VarWalker algorithm [46] to combine PPI network and mutation data identified by next-generation sequencing (NGS) to build consensus networks for identifying cancer driver genes. While VarWalker was proposed for predicting driver mutations, it provides a general framework for prioritizing target genes from high-throughput sequencing data assisted by PPI network. VarWalker evaluates candidate target genes (i. e. , mutation gene in cancer or DmMG in this work) by assessing their topological significance using random networks which hold the same degree distribution with the PPI network. However, utilizing only the topological characteristics may remove functionally significant candidate target genes. Also, the filtering result is not steady because it will remove different candidate genes for the same target gene when using different random networks. That is, VarWalker is not robust enough. We propose in this paper an improved strategy to evaluate both the topological and functional significance of candidate DmMGs in a more robust and efficient way, and the approach is detailed in the Materials and Methods section. To compare the robustness and efficiency of m6A-Driver and VarWalker in filtering candidate genes, we applied the two methods on 100 genes randomly selected from the DmMGs in KD-METTL3 dataset to filter their candidate genes using two different sets of random networks. Each set contains 100 networks which hold the same topological property of PPI network. A more robust algorithm should remove a consistent set of genes in two random network sets. As is shown in Table 1, m6A-Driver only removed 1 different candidate gene when using different random network sets, whereas VarWalker removed 40 different candidate genes. This result demonstrates that m6A-Driver is more robust in filtering candidate genes. It is not surprising to also notice that some of the removed genes by VarWalker have significant biological functional connections with the seed (the DmMG) in the PPI network. Moreover, m6A can filter candidate genes in a more efficient way. VarWalker needs to perform RWR for each DmMG in each of the 100 random networks to compute the reoccurrence frequency of the candidate genes for calculating the p-value. In contract, m6A-Driver only needs the degree of a candidate gene and the degree of the seed gene to calculate its empirical p-value. To validate m6A-Driver, we applied it to the four different case-control MeRIP-seq datasets: KD-FTO, KD-METTL3, KD-METTL14 and KD-WTAP. KD-FTO includes 9 RSs, based on which an FTO knockdown mDrNet (S2 Fig) was built. The network consists of 1,832 mDrGenes and 21,506 edges, with the maximal connected sub-graph containing 1,787 mDrGenes, implying that there exist dense interactions among mDrGenes. KD-METTL3, KD-METTL14 and KD-WTAP all include 4 RSs, based on which the corresponding context-specific mDrNets (S3–S5 Figs) were constructed by m6A-Driver. KD-METTL3 mDrNet contains 1,352 mDrGenes and 8,235 edges, with the maximal connected sub-graph including 1,339 mDrGenes; KD-METTL14 mDrNet consists of 1,251 mDrGenes and 8,452 edges, with itself being the maximal connected sub-graph; KD-WTAP mDrNet has 375 mDrGenes and 1,980 edges, which is also its maximal connected sub-graph. Similar to KD-FTO network, most mDrGenes in each of the 3 networks are interacting with each other very closely, implying again that the predicted mDrGenes have highly relevant functions. We next examined the characteristics of the predicted mDrGenes. We first investigated the differential methylation of the mDrGenes. An mDrGenes is defined as a hyper mDrGenes if its most differentially methylated site is hyper-methylated, but otherwise defined as a hypo mDrGenes if its most differentially methylated site is hypo-methylated. We counted the number of hyper and hypo mDrGenes (Fig 2). As expected, mDrGenes in KD-FTO are mostly hyper-methylated, whereas those in three other methylase knockdown datasets are more hypo-methylated. This result is consistent with the fact that FTO is an m6A demethylase, but METTL3, METTL14, and WTP are elements of m6A methyltransferase complex. We also calculated the average number of DmM sites for per gene and found that on average, an mDrGene harbors more than one DmM sites (Table 2). It is interesting that mDrGenes in KD-FTO harbor more DmM sites than the other 3 datasets and KD-METTL3 mDrGenes harbor the least number of DmM sites on average. We then investigated the DmM site distribution using the Guitar R/Bioconductor package [47] in an mDrGene transcript (Fig 3). Overall, the distributions for the 4 datasets are very similar, where DmM sites are mostly enriched around the stop codon and are distributed more in 3' UTR and CDS, which is consistent with the reported results in the literature [6,7]. Furthermore, we obtained the sequence motifs of DmM sites in mDrGenes for each of the four datasets using MEME-ChIP webserver [48] (Fig 4). The reported RRACH m6A motifs [6,7] was top ranked in KD-FTO and KD-METTL3, whereas the most enriched motifs in KD-METTL14 and KD-WTAP are similar to the binding motifs of SRSF1 and SRSF9. Interestingly, SRSF1 and SRSF9 are components of the SRSF protein that is involved in splice site selection in alternative splicing [49]. We asked if mDrGenes are more functional relevant. To test this, we examined the functionally significance between mDrGenes and DmMGs predicted by exomePeak. We performed GO [50] enrichment analysis using DAVID (Database of Annotation, Visualization and Integrated Discovery) [51] and then compared the enrichment degrees of the top enriched biological processes (BP, Fig 5). Since a larger testing gene set tends to lead to a smaller enriched p-value when performing DAVID, to make the comparison fair, we balanced the scale of mDrGenes and DmMGs before enrichment analysis. For KD-FTO dataset, the scale of mDrGenes is larger (Fig 6), so we randomly removed some mDrGenes to make the scales the same and then performed the enrichment analyses for 10 times to calculate an average pBonferroni for each enriched term. After also performing enrichment analysis on DmMGs, we compared the pBonferroni of top 20 enriched terms for DmMGs and scaled mDrGenes. Since, DmMGs have a larger scale for KD-METTL3, KD-METTL14 and KD-WTAP datasets, to balance the scale, we selected DmMGs that harbor top differently methylated DmM sites. Then, we performed enrichment analyses and compared the pBonferroni of top 20 enriched terms for the mDrGenes and scaled DmMGs. The result shows that mDrGenes are more significantly enriched than DmMGs in all the top enriched biological processes, demonstrating that mDrGenes are more functional relevant than DmMGs. To further investigate the biological significance of mDrGenes, we evaluated the differential expression (DE) of mDrGenes and DmMGs. Conceivably, a gene set is likely to be more functionally important if it has more differentially expressed genes (DEGs) and/or its DEGs are more differentially expressed. To this end, we applied DESeq2 [52] to the input replicates of treated and untreated samples and determined a gene to be DEG if the adjusted p-value is less than 0. 05. We first examined the percentage of DEGs in mDrGenes and DmMGs in the four datasets. We found that there are very few DEGs in mDrGenes and DmMGs for both KD-FTO and KD-METTL14 dataset, and thus the percentages of DEGs in mDrGenes and DmMGs are very low for these two datasets. Not surprisingly, no significant differences between the percentages of DEGs in mDrGenes and DmMGs can be observed (Fisher’s test, see Table 3 for details). In contrast, much more mDrGenes and DmMGs are differential expressed in KD-METTL3 and KD-WTAP, and the percentages in mDrGenes are significantly higher than those in DmMGs (Fisher’s test, see Table 3 for details). We next compared the degree of DE, which is represented by the negative log10 (FDR) calculated by DESeq2 (Fig 7). The result shows that the DE degrees of mDrGenes are also higher than those of DmMGs. The only exception is the FTO KD experiment, in which there are nearly no differential expression genes. Taken together, we can conclude that mDrGenes are likely to include more DEGs than DmMGs and their degree of DEs are also likely to be higher. Functional enrichment analyses were carried out on the 4 mDrNets to help reveal the biological processes regulated by the 4 enzymes at epitranscriptomic layer of gene regulation. The results obtained using DAVID reveal a significant enrichment of multiple m6A-related pathways annotated by either the Kyoto Encyclopedia of Genes and Genomes (KEGG) [53] or Gene Ontology (GO) biological process (BP) domains (Figs 8–11). For KD-FTO, since FTO is a demethylation enzyme, we expected to observe mainly hyper-methylation. However, m6A-Driver did report several hypo-methylation mDrGenes, suggesting a potentially direct or indirect mode of FTO regulation to also enhance m6A. We next examined the functional relationship between the hyper- and hypo-methylated genes and found that there was little overlapping between their enriched functions (Fig 8 and Fig 9, see S6 Fig and S7 Fig for detail). To keep consistency with paper [27], in which the KD-FTO data is published, we adopt the whole reference genome as the control data set of enrichment analysis. The hyper-methylated mDrGenes are clearly linked closely to neurological processes and neuro signaling pathways. Several significantly enriched terms annotated by GO BP are synapse and neuron signaling transmission (132 hyper mDrGenes, pBonferroni = 1. 79×10−18), synaptic transmission (43 hyper mDrGenes, pBonferroni = 1. 77×10−13), transmission of nerve impulse (50 genes in hyper mDrGenes, pBonferroni = 5. 58×10−14). They are also likely associated with neuron differentiation (65 hyper mDrGenes, pBonferroni = 1. 81×10−11) and neuron development (53 hyper mDrGenes, pBonferroni = 3. 00×10−11) as well as embryonic development (54 hyper mDrGenes in utero embryonic development, pBonferroni = 2. 31×10−13 and 67 hyper mDrGenes in chordate embryonic development, pBonferroni = 2. 42×10−11), which may be another proof of RNA methylation involved in steering stem cell pluripotency [34–36]. In contrast, the hypo-methylated mDrGenes are more related to metabolic processes (79 hypo mDrGenes in protein catabolic process, pBonferroni = 2. 41×10−18 and 85 hypo mDrGenes in macromolecule catabolic process, pBonferroni = 1. 88×10−17) and cell cycle (65 hypo mDrGenes, pBonferroni = 9. 29×10−12). In addition, the hypo mDrGenes are enriched in Spliceosome (17 hypo mDrGenes, pBonferroni = 8. 25×10−4), which is also a KEGG term enriched in KD-METTL3, KD-METTL14 and KD-WTAP data (Fig 11), implicating a potential role of m6A in mRNA splicing. Note that WTAP itself is also splicing factor. However, this result suggests that WTAP might also regulate splicing in an m6A dependent fashion. Taken together, our predicted mDrGenes confirm the demethylation role of FTO but may suggest a direct or indirect role of FTO in promoting m6A. Functional enrichment suggests that these two modes of FTO function are involved in distinct biological processes and pathways. Another interesting finding is that both hyper and hypo mDrGenes are enriched in cancer related pathways including Chronic myeloid leukemia (20 hyper mDrGenes, pBonferroni = 3. 27×10−6; 12 hypo mDrGenes, pBonferroni = 2. 12×10−3) and Glioma (16 hyper mDrGenes, pBonferroni = 8. 01×10−5; 12 hypo mDrGenes, pBonferroni = 4. 86×10−4) (Fig 9, see S7 Fig for detailed information). As FTO-KD data is extracted from mouse brain, so we also do enrichment analysis using the brain tissue specific expressed genes as control data to check whether the pathways we find above are really influenced by m6A methylation or by tissue specific expression. Brain tissue specific expressed genes are defined here as genes who have a RPKM value over 1 in at least half of the input samples, including treated and untreated ones. As is shown in S8 and S9 Figs, the results are similar to make the whole reference genome as control data set, but the enriched pathways get a bigger q-value due to the size reducing of control datasets, such as synaptic transmission (43 hyper mDrGenes, pBonferroni = 1. 09×10−5), transmission of nerve impulse (50 genes in hyper mDrGenes, pBonferroni = 1. 98×10−5), neuron differentiation (65 hyper mDrGenes, pBonferroni = 1. 89×10−8) and Pathways in cancer (64 hyper mDrGenes, pBonferroni = 5. 12×10−6; 40 hypo mDrGenes, pBonferroni = 3. 51×10−2). These results show that the mDrGenes enriched pathways are really influenced by m6A methylation. To further study the dynamics of m6A methylation, we applied the enrichment analysis to mDrGenes predicted in KD-METTL3, KD-METTL14 and KD-WTAP and compared the similarity and difference of their enriched GO biological processes (Fig 10, and also see S10 Fig for more details) and KEGG pathways (Fig 11, and also see S11 Fig for more information). In this case, we chose to perform enrichment on all mDrGenes in these three datasets instead of analyzing hyper- and hypo-methylated mDrGenes separately, because there are little hyper mDrGenes in these 3 datasets. There are significant overlapping biological processes among these 3 sets of mDrGenes, but also exist enzyme specific functions (Fig 10). The common biological processes include cell cycle (143 mDrGenes in KD-METTL3, pBonferroni = 2. 52×10−20,123 mDrGenes in KD-METTL14, pBonferroni = 6. 71×10−15,51 mDrGenes in KD-WTAP, pBonferroni = 2. 35×10−10), regulation of transcription (324 mDrGenes in KD-METTL3, pBonferroni = 7. 76×10−16,315 mDrGenes in KD-METTL14, pBonferroni = 2. 37×10−19,103 mDrGenes in KD-WTAP, pBonferroni = 1. 88×10−7) and positive regulation of molecular metabolism, e. g. , positive regulation of macromolecule metabolic process (137 mDrGenes in KD-METTL3, pBonferroni = 5. 19×10−14,126 mDrGenes in KD-METTL14, pBonferroni = 8. 12×10−13,55 mDrGenes in KD-WTAP, pBonferroni = 9. 27×10−11). What is interesting is that the overlapping functions between METTL3 associated mDrGenes and WTAP associated mDrGenes are positive regulations of metabolism and gene expression. In contrast, the overlapping functions between METTL3 targeted mDrGenes genes and METTL14 mDrGenes are mainly negative regulation of metabolism and gene expression. The METTL14 mDrGenes are also strongly enriched in splicing, especially, e. g. , RNA splicing (76 mDrGenes, pBonferroni = 2. 71×10−22), RNA splicing, via transesterification reactions (43 mDrGenes, pBonferroni = 1. 38×10−13) and RNA splicing, via transesterification reactions with bulged adenosine as nucleophile (43 mDrGenes, pBonferroni = 1. 38×10−13). It is consistent with that the most enriched motif in KD-METTL14 is similar to the binding motifs of SRSF1 and SRSF9, two factors involved in alternative splicing. These suggest a potential role of METTL14 in regulating splicing via m6A. In contrast, WTAP mDrGenes are enriched specifically in chromatin modification (26 genes, pBonferroni = 1. 46×10−8), whereas METTL3 mDrGenes may influence the development of protein complex, e. g. , protein complex assembly (96 mDrGenes, pBonferroni = 1. 86×10−14) and protein complex biogenesis (96 mDrGenes, pBonferroni = 1. 86×10−14). Comparison of top KEGG pathways enriched in the three mDrGenes sets also revealed common as well as methylase specific functions (see Fig 11). Particularly, the consensus functions include cell cycle (38 mDrGenes in KD-METTL3, pBonferroni = 2. 19×10−11,27 mDrGenes in KD-METTL14, pBonferroni = 8. 64×10−7,16 mDrGenes in KD-WTAP, pBonferroni = 1. 08×10−5), spliceosome (34 mDrGenes in KD-METTL3, pBonferroni = 8. 9×10−9,26 mDrGenes in KD-METTL14, pBonferroni = 3. 51×10−6,14 mDrGenes in KD-WTAP, pBonferroni = 2. 06×10−4) and pathway in cancer (66 mDrGenes in KD-METTL3, pBonferroni = 2. 68×10−10,57 mDrGenes in KD-METTL14, pBonferroni = 8. 58×10−10,30 mDrGenes in KD-WTAP, pBonferroni = 5. 46×10−7) especially Chronic myeloid leukemia (28 mDrGenes in KD-METTL3, pBonferroni = 8. 66×10−11,21 mDrGenes in KD-METTL14, pBonferroni = 2. 52×10−7,12 mDrGenes in KD-WTAP, pBonferroni = 2. 58×10−5). The significant overlapping pathways between METTL3 and METTL14 include Glioma (20 mDrGenes in KD-METTL3, pBonferroni = 1. 37×10−6,13 mDrGenes in KD-METTL14, pBonferroni = 1. 85×10−3), suggesting that these mDrGenes may be used as biomarkers of glioma. We also notice that METTL3 mDrGenes are specifically enriched in Melanoma (18 mDrGenes, pBonferroni = 1. 32×10−4). A recent study have demonstrated that mutations within intron 8 of FTO leads to increased melanoma risk [29], suggesting a link between m6A and melanoma. To help reveal the pathways potentially relevant to different modes of m6A functions, we checked the overlaps between the enriched pathways in hyper and hypo mDrGenes in KD-FTO and mDrGenes in KD-METTL3, KD-METTL14 and KD-WTAP (S12 and S13 Figs). All 5 groups of mDrGenes are enriched in cell cycle and pathways in cancer, including especially Chronic myeloid leukemia. This further suggests that m6A is related to cancer. The overlapping pathways between the hyper-mDrGenes in KD-FTO and those in METTL3/METTL14 are mainly related to transcription including regulation of transcription and regulation of gene expression. Indeed, it has been shown that m6A recruits YTHDF2 protein to regulate mRNA stability [54]. In contrast, the overlapping pathways between the hypo-mDrGenes in KD-FTO and those in METTL3/METTL14 are related RNA splicing. Interestingly, nuclear m6A-binding protein YTHDC1 is shown to promote exon inclusion of targeted mRNAs through facilitating mRNA binding of splicing factor SRSF3 [55,56]. We further examined the mDrNets and their subnetworks associated with the enriched biological processes. Several sub-mDrNets for KD-FTO including intracellular signaling cascade, synaptic transmission BP category, and transmission of nerve impulse, are shown in Fig 12. They are consistent with our hypothesis that mDrGenes in the same biological process are interacted with each other closely. That is also the reason why m6A-Driver can identify mDrGenes that might not be consistently identified in most RSs due to the biological variance, but have significant biological functions. This underscores the advantage of m6A-Driver in addressing variance among different replicates for predicting mDrGenes. For KD-METTL3, KD-METTL14 and KD-WTAP, since they form the m6A methyltransferase complex, we integrated the 3 mDrNets and examined the subnetworks associated with the enriched pathways (Fig 13). Similar to KD-FTO, the mDrGenes enriched in the same pathway are closely connected and many mDrGenes undetected as differential m6A genes in all RSs are also identified. What is also interesting to notice is that the enriched pathways common in three datasets could be resulted from very different mDrGenes for each dataset, suggesting the 3 m6A methylases collaboratively regulate the same pathway through different mechanisms. Because many m6A sites are detected in 3’UTRs that also contain microRNA binding sites. It will be helpful to further examine if the mDrGenes are also enriched in certain microRNA with significant functions. To test this, we download the microRNA-target information from miRTarBase, a database curates experimentally validated microRNA-target interactions [57]. Then we performed Fisher’s test to test whether the m6A driven genes are enriched in targets of certain microRNA families. Interestingly, although most mDrGenes are targeted by microRNAs (60% in KD-FTO, 96% in KD-METTL3,97% in KD-METTL14 and 97% in KD-WTAP), not many microRNAs have targets enriched in mDrGenes. For KD-FTO, there are only 2 microRNAs have p-value < 0. 05, and for KD-METTL3, KD-METTL14 and KD-WTAP, there are only 1,1, 0, separately (Table 2). The information of all targeted mDrGenes by microRNA is included in supplementary material (S1–S4 Texts).
The MeRIP-seq technology significantly advances the study of m6A, enabling profiling m6A methytranscriptome for specific cell conditions. However, existing algorithms focus mostly on predicting m6A sites from MeRIP-seq data. Although they are powerful tools for MeRIP-seq data analysis, they cannot directly assess the functional importance of these sites and associated genes. To address this shortcoming, we proposed in this paper m6A-Driver, a novel algorithm for detecting m6A-driven genes and their interaction network. m6A-Driver utilizes protein-protein interaction networks to identify functional meaningful differentially m6A methylated genes and overcomes the biases in predicting functional enrichment of sites due to different sample replicates. The comparison on the p-values of top enriched biological processes in the prediction results of m6A-Driver and exomePeak demonstrates that m6A-Driver could identify mDrGenes that are more functional relevant. In terms of the algorithm, comparison with VarWalker, an algorithm for predicting mutation driver genes, shows that m6A-Driver is computationally more efficient and can produce topological and biologically more robust predictions. Furthermore, m6A-Driver generates a condition-specific m6A-driven network that reveals the detailed functional circuitry underlying the biological condition. The results on the FTO, METTL3, METTL14 and WTAP knockdown data demonstrated that m6A-Driver can address the sample bias in MeRIP-Seq data and identify functional relevant mDrGenes in a robust and efficient fashion. m6A-Driver predicted several significant biological progresses and pathways associated with each knockdown dataset and constructed four mDrNet separately regulated by the four m6A (de) methylases. Functional enrichment analysis across the four networks showed the involvement of m6A in a diverse processes and pathways including regulating cancer, transcription, and slicing, many of which have been reported in the literature. The presented results in this paper demonstrate that m6A-Driver is an effective and reliable approach to identify functionally relevant m6A-driven genes and networks from MeRIP-Seq data. We want to point out that the inherent technical limitations of MeRIP-Seq can lead to increased false positives in the predicted mDrGenes and their functional networks. As pointed out earlier in the paper, the sample bias and the impact of library “size factor” from different conditions can negatively impact the quality of the data. How to normalize MeRIP-Seq samples from different conditions is still an open topic that requires additional research. Many useful normalization methods such as the" geometric mean" approach proposed in DESeq can provide valuable guidance. Nevertheless, the results from this and other papers showed that with careful processing and modeling of the current MeRIP-Seq data, many important functions of m6A can be predicted. With the continuing improvement of the MeRIP-Seq and related technologies, we expect that m6A-Driver should produce increasingly more accurate predictions. In the current m6A-Driver algorithm, identification of mDrGenes relies on a reference network and a threshold to determine the candidate genes. So far, only PPI network is considered, and no reference networks for noncoding RNAs are included. As m6A is also prevalent in long noncoding RNAs (lncRNAs), constructing a noncoding RNA interaction network and extending m6A-Driver algorithm to lncRNAs will enable the study of m6A associated functions in lncRNAs. As a future work, we can use the score in step 4 to help construct the consensus network. For a DmMGs network of a biological replicate sets built by step 3, we can use peak calling score to transfer the network to a weighted network. The weight of an edge in the network could be calculated as w=−lg (x) −lg (y) 2, where x and y denote the peak calling FDR (adjusted p-value) of the two gene nodes. Then all edges of the DmMGs networks are pooled together. For each edge, we can calculate the consensus score as s=∑i=1mwi, where m is the number of replicate sets and wi is the weight of the edge in network i. The final consensus m6A driven gene network can be built by setting a threshold on s according to its distribution or the scale of the consensus network we need.
FTO knockdown (KD-FTO) dataset is a MeRIP-Seq data [27] from the wild-type littermate as well as FTO knockdown mice samples. There are 9 sets of biological replicates in this cohort, and each biological replicate set (RS) contains two IP samples from a FTO knockdown mouse and the wild-type littermate, respectively, as well as two corresponding input samples from the two mice samples. In the original work [27], Hess et al. demonstrated that FTO-knockdown mice have impaired dopamine release, reduced dopaminergic receptor responses and an altered locomotor response to cocaine, which are related to specific m6A mRNA demethylation regulated by FTO. METTL3, METTL14 and WTAP knockdown MeRIP-Seq datasets are from a recent study which reveals that m6A regulates mRNA stability [44]. Each of the three cohorts consists of 4 RSs, each of which contains two IP samples from gene knockdown (treated) HeLa cells and untreated HeLa cells, respectively, and two corresponding input samples. The study revealed that knockdown of METTL3 led to decrease the binding of YTHDF2 to its targets, and increase the stability of its target RNAs similar to that of YTHDF2 knockdown. The reference network, PPI network, is built from the most recent version of PPI data from BioGRID (release 3. 4. 128, compiled on August 25th, 2015) [45]. Based on the binary interactions, we removed the isolated and self-interactions proteins to establish a PPI network with a total of 16,062 proteins and 152,676 interactions. m6A-Driver predicts mDrGenes and mDrNet from MeRIP-Seq data, where mDrGene is a gene whose mRNA harbors at least one DmM site in a biological context of interest and whose function is relevant to the context. m6A-Driver first divides the MeRIP-Seq data into several RSs, each containing 2 paired samples, an IP sample paired with an input sample under the treated condition and another pair under the untreated condition. Then, m6A-Driver predicts mDrGenes and mDrNet by the following four steps. The workflow is shown in Fig 1.
Powered by methylated RNA immunoprecipitation sequencing (MeRIP-Seq) technology, recent studies have revealed a new mode of post transcriptional regulation mediated by mRNA N6-methyladenosine (m6A). Currently, the analysis of m6A focuses mostly on prediction of m6A sites as well as differential m6A methylation, and systematic approach for predicting m6A functions is yet to emerge. We develop here m6A-Driver, the first network-based approach, to identify m6A-driven genes and their associated networks, whose functional interactions are likely to be actively modulated by m6A methylation under a specific condition. Our test results showed that m6A-Driver can robustly and efficiently identify m6A-driven genes that are functionally more enriched and associated with higher degree of differential expression than differential m6A methylated genes. m6A-Driver is an effective and reliable approach to identify functionally relevant m6A-driven genes and networks from MeRIP-Seq data.
genetic networks
gene regulation
applied mathematics
protein interaction networks
enzymes
enzymology
simulation and modeling
algorithms
methylation
mathematics
network analysis
enzyme chemistry
research and analysis methods
computer and information sciences
proteins
enzyme regulation
gene expression
chemistry
proteomics
methyltransferases
biochemistry
gene identification and analysis
genetics
biology and life sciences
chemical reactions
physical sciences
m6A-Driver: Identifying Context-Specific mRNA m6A Methylation-Driven Gene Interaction Networks