Limits in the detection of m6A changes using MeRIP/m6A-seq

Many cellular mRNAs contain the modified base m6A, and recent studies have suggested that various stimuli can lead to changes in m6A. The most common method to map m6A and to predict changes in m6A between conditions is methylated RNA immunoprecipitation sequencing (MeRIP-seq), through which methylated regions are detected as peaks in transcript coverage from immunoprecipitated RNA relative to input RNA. Here, we generated replicate controls and reanalyzed published MeRIP-seq data to estimate reproducibility across experiments. We found that m6A peak overlap in mRNAs varies from ~30 to 60% between studies, even in the same cell type. We then assessed statistical methods to detect changes in m6A peaks as distinct from changes in gene expression. However, from these published data sets, we detected few changes under most conditions and were unable to detect consistent changes across studies of similar stimuli. Overall, our work identifies limits to MeRIP-seq reproducibility in the detection both of peaks and of peak changes and proposes improved approaches for analysis of peak changes.

Intersection Size Replicates  (the number of DRAC motifs divided by peak length, with mean indicated above the violin plots) across methods.
Boxes within the violin plots span the 1 st to 3 rd quartiles, with medians indicated. Whiskers show the minimum and maximum points within ±1.5x the interquartile distance from the boxes, with points beyond that range shown as outliers. b) MeRIP-RT-qPCR analysis of a subset of peaks called by single tools to test for METTL3-METTL14 dependence. Huh7 cells were treated with siRNAs targeting METTL3 and METTL14 (siMETTL) or non-targeting control siRNA (siCTRL) and MeRIP-RT-qPCR was performed using primers under uniquely called peaks. HPRT1 (known unmethylated transcript) and RNU6 (known to be methylated by METTL16) were negative controls.
Positive controls were previously published known methylated transcripts. c) The percent of genes with ≥ 1 m6A peak for various mean coverage bins (0-0.05, 0.05-0.1, etc.). Results are summarized across the experiments included in Figures 1 and 2. Boxplot center lines show medians and whiskers 1.5x interquartile range. d) Histograms of peaks detected at increasing median input levels across replicates confirms that peak detection requires adequate input coverage (here, few peaks are called at input levels below ten). e) METTL3/14dependence of a random selection of peaks detected in low expression < 10 reads vs. higher expression > 10 reads sites from Huh7 data. MeRIP-RT-qPCR analysis of relative m 6 A level RNA harvested from Huh7 cells that were treated with siCTRL or siMETTL3/14. Values are the mean ± SD of 3 replicates. * p < 0.05, ** p < 0.01, *** p < 0.001 by unpaired Student's t test. n.s = not significant. f) The correlation between log2 fold IP/input read counts between two random replicates. g) The enrichment of DRAC motifs calculated as # DRAC motifs/peak length for peaks detected by X replicates shows little change for peaks detected in multiple replicates, with the number of peaks detected by X replicates shown below. The lines in the violin plots indicate the median. a,d,f,g show data from mouse cortex (left) and Huh7 cells (right) at baseline conditions.

Supplementary Figure 2. Results using the MeTDiff peak caller. a) Peak detection using MeTDiff between
studies that used the same cell type shows variable overlap. Overlap was calculated as the percent of peaks detected in Experiment 1 with an overlap of ≥ 1 base pair with peaks from Experiment 2 (compare to Figure 2a), b) Peak detection across tissue and cell types shows samples from the same study cluster better together than samples from the same tissue (compare to Figure 2b). Median overlap was 55%.  Table 1) to detect differential methylation in negative controls between two groups at baseline conditions and positive controls in which methylation processes were disrupted with respect to baseline conditions (Supplementary Table 3 Table 4). The number of peaks detected as changed in the original published analyses are compared to the number of peaks with FDR-adjusted p-values < 0.05 in our reanalysis using DESeq2, edgeR, or QNB, and taking the union of results from these three tools with additional filters for log2 fold difference in peak and gene changes of ≥1 and peak read counts ≥10 across all replicates and conditions ("filtered") (compare to Figure 4a), b) The same, without a threshold for peak read counts (compare to Supplementary Figure 7b).   Dashed lines indicate differences in log2 fold change of -1, 0, and 1 and medians are summarized below the distributions. Boxes within the violin plots span the 1 st to 3 rd quartiles, with medians indicated. Whiskers show the minimum and maximum points within ±1.5 times the interquartile distance from the boxes, with points beyond that range shown as outliers.

Supplementary Figure 7. a)
Coverage for a gene involved in pluripotency that we detected as significantly differentially methylated with activin-NODAL inhibition by SB431542 (SB). The peak did not pass the filter for minimum input reads. b) Detected m 6 A(m) changes in thirteen published data sets that measured m 6 A(m) peak changes between two conditions (Supplementary Table 4). The number of peaks detected as changed in the original published analyses are compared to the number of peaks with FDR-adjusted p-values < 0.05 in our reanalysis using DESeq2, edgeR, or QNB, and taking the union of results from these three tools with additional filters for log2 fold difference in peak and gene changes of ≥1 ("filtered"), without an additional filter for input read counts. c) The sole peak detected as changed in the Huang et al. (2019)  data in which we did not detect any significant changes in methylation, with or without a filter for input reads. In general, low input coverage in the first two replicates from this study appeared to contribute to the elevated standard deviations in coverage. f) The correlation between changes in enrichment (Δ peak -gene log2 fold change between conditions) at the same sites between experiments that studied similar conditions. a,c-e) Lines show the mean coverage across three replicates for input and IP samples, while shading shows the standard deviation. Peaks detected as significantly changed are highlighted in yellow and coding sequences are shown in grey.