m6A-seq data analysis

m⁶A peak identification was performed as previously described (32). In brief, we trimmed the adaptor and low quality reads using Cutadapt (33) (-e 0.3 –minimum-length 25; –trim-n -q 20,20) and fastx toolkit (fastx_trimmer -f 6; -t 5 -m 20). rRNA reads were then removed using SortMeRNA (34). Next, clean reads were mapped to the human genome (hg19) using TopHat2 (version 2.1.0) (35). Enriched peaks were identified by scanning each gene using 100-nt sliding windows, and calculating an enrichment score for each sliding window (winscore).

MeanWinIP and MeanWinControl are the mean coverage for each window for immunoprecipitation and input control, respectively. MedianGeneIP and MedianGeneControl are gene median coverages for immunoprecipitation and input control, respectively. Windows with RPKM ≥ 10 in the IP sample, enrichment score ≥2 in genes with RPKM in the input sample ≥1 were defined as enriched windows. Last, only the peaks with winscore ≥2 in at least two samples of a tissue type were considered as real m⁶A peaks.

To generate the metagene profile of m⁶A site distribution across transcripts, we first determined the number of bins that need to be divided for a given gene based on relative lengths between 5′UTR, CDS and 3′UTR of the human transcriptome (GENCODE v26). The relative lengths between 5′UTR, CDS and 3′UTR are 10:50:40, thus for each gene, 10, 50 and 40 bins of equal length were made for 5′UTR, CDS and 3′UTR, respectively. Next, for each m⁶A site, we assigned it to the longest isoform of the corresponding gene and determined which bin it is located in. Last, the number of sites for each bin was summarized and the curve was fitted with polyfit.

The mapping statistics of all datasets were summarized in Supplementary Table S1.