m6A peak calling, differential peak calling, and motif analysis

To identify regions in which m⁶A modifications occur, we used the peak calling algorithm MACS2 (version 2.1.0.20140616) [62] on m⁶A-seq, using the corresponding RNA-seq as input (background). Studies have shown that MACS2 is a suitable tool to analyze MeRIP-seq data. For example, work by Liu et al. [63] compared MACS2 and exomePeak and showed a significant correlation between the results obtained on MeRIP-seq datasets using the two techniques. And a more recent study by McIntyre et al. [64] drew the same conclusions, stating that most people use MACS2 to analyze MeRIP-seq and they found it to be a reliable tool. In fact, these authors found MACS2 outperforms exomePeak. Another study by Anatanaviciute et al. [65] also shows that MACS2 was a better performer than exomePeak. The MACS2 callpeak function was run with the following parameters: –nomodel,–extsize 150, -p 5e-2, and –g mm. Therefore, candidate m⁶A peaks were identified as an enrichment of reads upon pull down with the m⁶A antibody compared to background. To identify differential peaks, the MACS2 differential binding events program (bdgdiff) with parameters -g 20 and -l 120 was employed. The consensus sequence motifs enriched in m⁶A peaks were identified by using MEME (version 5.1.1) [66]. The integrative genomics viewer (IGV) tool was used for visualization of m⁶A peaks along the whole transcript [67].