m6A peak calling and motif analysis

RNA m⁶A-modified regions, also called m⁶A peaks, were identified using exomePeak software (version 2.7.0) with FDR (false discovery rate) < 0.05 [66, 67]. The consensus sequence motifs enriched in m⁶A peaks were identified by HOMER [68]. For further comparison of m⁶A modification between samples, we used coverageBed of BEDToods (version 2.26.0) [69] with the “-s”, “-splited”, “-counts”, and “-F 0.50” parameters to calculate the read count of each peak. The “IP FPKM”, “input FPKM”, and “Enrichment score” of peaks were calculated as following methods:

where a_i,j denotes the peak “IP FPKM” value of the ith methylation peak in the IP sample from the jth biological sample; A_i,j denotes the total reads mapped to the ith methylation peak in the IP sample from the jth biological sample; B_j denotes the total unique reads mapped to the mouse reference (mm10) in the IP sample from the jth biological sample; C_i,j denotes the length (base) of the ith methylation peak in the IP sample from the jth biological sample; b_i,j denotes the peak “input FPKM” value of the ith methylation peak in the input sample from the jth biological sample; D_i,j denotes the total reads mapped to the ith methylation peak in the input sample from the jth biological sample; E_j denotes the total unique reads mapped to the mouse reference (mm10) in the input sample from the jth biological sample; c_i,j denotes the peak “enrichment score” value of the ith methylation peak in the IP/input sample from the jth biological sample.

M⁶A peaks that satisfied 1) exomePeak FDR value < 0.05, 2) m⁶A peak “IP FPKM” value > 1, and 3) m⁶A peak “enrichment score” value > 1.5 were used for further comparative analysis. The continuously methylated RNAs (CMRs) were defined as RNAs containing at least one m⁶A peak in all the four samples, while the temporal-specific methylated RNAs (SMRs) were the RNAs with m⁶A modification only in one sample, and m⁶A peaks of SMRs were defined as specific m⁶A peaks.

By merging all the m⁶A peaks of four developmental stages using mergeBed of BEDTools with “-s”, “-c”, and “-o” parameters, we identified “ON” or “OFF” RNA m⁶A switches during cerebellar development. All m⁶A peaks absent in the former stage but present in the later stage were called “ON” m⁶A switches, while all m⁶A peaks displaying changes in the opposite direction were called “OFF” m⁶A switches. In addition, the differentially methylated RNAs (DMRs) between wild-type and Alkbh5 knockout mice were further identified using exomePeak software with default parameters, and high-confidence DMR m⁶A peaks matching the criteria of FDR < 0.05 and log2 (fold change) > 1 or < − 1 were selected for further analysis.

For validation analysis of the distribution patterns of stage-specific m⁶A peaks which were obtained by using the exomePeak software, peak calling analysis was performed in parallel using the MACS2 software for the P7 and P60 samples with the default parameters [54, 70].

The integrative genomics viewer (IGV) tool was used for visualization of m⁶A peaks along the whole transcript [71, 72]. The heatmaps of the enrichment score for m⁶A peaks were plotted using the pheatmap package in R [73].