We assess the effect of experimental design on the power of DMR detection purely based on simulations, where the whole procedure is divided into two components. First, magpie preprocesses .bam files from MeRIP-seq sequencing and obtains read counts in candidate regions from all samples (Fig. 1), where candidates are identified with conditional binomial tests. With the counts from the identified candidate regions, magpie simulates count matrices for both IP and Input samples with a Gamma-Poisson model. Parameters involved are estimated from the candidates to mimic the actual MeRIP-seq data in aspects of marginal distribution read counts, and the distribution of biological dispersion in methylation levels (Fig. 1). With data simulated, we then evaluate power and error rates on them (Fig. 1). The two components, Gamma-Poisson simulation and power assessment, are independent so that magpie allows the assessment on data by different simulation strategies.

Here we describe how magpie simulates the MeRIPseq count matrices given existing real MeRIP-seq data from different conditions. magpie processes .bam files by splitting the transcriptome into bins, aggregating read counts, and testing for significance of IP enrichment over Input. Using a bump-finding algorithm, significant bins are combined into candidate regions. magpie then focus on these candidates in simulations, as other regions lack IP enrichment and biological relevance. Suppose there are in total N pairs of IP and Input samples from all conditions, and M candidate DMRs generated after preprocessing. Let Xij and Yij denote input and IP counts in candidate DMR i from sample j. We assume that Xij ∼ P oisson(sjxλixj ) and λixj ∼ Gamma(αixj , θi). Similarly, Yij ∼ P oisson(sjy λiyj ) and λiyj ∼ Gamma(αiyj , θi). Here sjx and sjy represent the normalizing factors for input and IP samples, such y as the library sizes. λixj and λij are normalized poisy son rates. αixj , αij , and θi are the shape and scale parameters of corresponding gamma distributions. Given y above assumptions, naturally λixjλ+ijλiyj ∼ Beta(αiyj , αixj ).

y αij y , φij = αixj +αij

1 αixj +αiyj +1 , then Further, denote μij = marginally,

Xij |μij , φij ∼ N B((1 − μij )(φi−j1 − 1), 1 +sjxsθjxiθi ), y sj θi ) Yij |μij , φij ∼ N B(μij (φi−j1 − 1), 1+sjyθi (1) In above equations, μij and φij represent the mean and dispersion of the methylation level for candidate region i in sample j.

We begin by simulating size factors, for which we directly use the values estimated from real data: sjx =

X xbj /median{X xb1, . . . , X xbN },

b b b sjy =

X ybj /median{X yb1, . . . , X ybN }

b b b where xbj and ybj are read counts in bin b from the jth Input and IP samples.

Next, for each candidate region i, magpie simulates a baseline methylation level μi or equivalently eαi in the structure of log μij = αi + zjT βi where zj contains the attributes of sample j, and βi represents corresponding effects. To do that, we randomly sample αi from one parametric distribution, or from its empirical distribution estimated from real data.

After simulating the baseline methylation, we simulate β0s for all regions. Because we can hardly know the i | 2 actual number of DMRs and their degree of differential methylation, specific settings are adopted based on both reasonable assumptions and empirical observations. First, magpie sets the proportion of DMRs as 10%, assuming that DM is present in only a small subset of regions in most experiments. Then, for nonDMRs, βi = 0. For DMR i, βi = βˆi if its estimated effect size βˆi is greater than the 50% quantile of all regions. Otherwise, βi ∼ U (1, 2). The exact quantile and the range of uniform distribution can be adjusted by users to explore a broad range of effect sizes.

For dispersion φij , we can simulate it again from a parametric distribution or sample from empirical distributions. To ensure the robustness, the empirical distribution can be estimated by TRESS from raw counts or by Beta-binomial regressions from normalized counts. Denote Y˜ij and Tij = X˜ij + Y˜ij as the normalized IP and total counts, the Beta-Binomial regressions are established as follows: ˜

Yij |Tij ∼ Bin(Tij , pij ) pij |μij , φi ∼ Beta(μij , φi), logit(μij ) = a˜i + zjT β˜i. (2) where μij and φi represent the mean and dispersion of methylation level. As noted, for the convenience of estimation, above Beta-binomial regressions (as well as TRESS) assume φij = φi for all j.

Empirically for the same real dataset, φ0 s estimated i by the Negative-binomial model in TRESS are usually greater than those estimated by Beta-binomial regressions. Without golden truth, our setting for φi relies on a data driven approach. Specifically, by comparing to the real data, magpie will calculate a KL-divergence for each of the synthetic counts by model in Eq. (1). Those φˆi resulting in a significantly lower KL divergence between simulated data and real data will be kept for final data generation. If there is no significant difference in KL-divergence, φˆi estimated by NB will be adopted. Lastly, we simulate the scale parameter θi in Eq. (1). Again, it can be simulated by parametric distributions, Simulated candidate count mplee matrix Sa siz . … .

… Sequednecpitnhg

After generating the simulated read counts in candidate DMRs, an existing software developed for MeRIP-seq is applied to detect DMRs. We implement an interface for calling TRESS and exomePeak2 (Tang et al., 2021) . Users can define other DMR detection methods and plug into the procedure. Each method reports test statistics, P-values and FDRs for all candidate regions. These results are then used for the downstream power assessment.

We adopted several evaluation metrics in the statistical power assessment for differential analysis using MeRIP-seq data. These metrics include classic criteria in hypothesis testings such as the false discovery rate (FDR), power, and precision. We also inspected the false discovery cost (FDC, defined below) and targeted power (Wu et al., 2015) , aiming to provide a comprehensive statistical power evaluation.

Because not all DMRs are of biological interest to us, especially those with low effect sizes, we introduce a cutoff Δ for the effect size β. Only those DMRs with |β| ≥ Δ are considered as ‘targeted DMRs’, which are of biological interest in research. We denote the number of non-DMRs, non-targeted DMRs, and targeted DMRs as R0, R1, and R2, respectively. Suppose Tr represents the testing result of region r, where Tr = 1 denotes the discovered DMR, and Tr = 0 otherwise. The confusion 3 5 3 5 7 10 matrix in the DMR detection is summarized in Table 1.

The false discovery rate (FDR) and precision are statistical metrics that jointly provide insights into the balance between true and false discoveries among the significant features. In practical experiments, such as MeRIPseq, researchers often focus on the top detected DMRs by their chosen methods for further analyses. As such, the proportion of true DMRs among the significants is essential to ensure the reliability of their downstream analyses. Thus, type I error control is often a prime task. In this context, FDR and precision are defined as E h PP0 i and E h P1P+P2 i, respectively. Statistical power is defined, naturally, as E h RP11++PR22 i. To investigate the power of detecting targeted DMRs that are biologically interesting with |β| ≥ Δ, targeted power is introduced and defined as E h RP22 i. To better illustrate the tradeoff between false positives and true positives, we propose an additional metric, False Discovery Cost (FDC), E h PP02 i, which is defined as the expected number of false positives per targeted true positive. The rationale behind this is straightforward: this cost is the expected number of false discoveries, per true discovery we are interested in.

Finally, our proposed evaluation framework allows for the examination of aforementioned metrics using simulations under various combinations of sample size, sequencing depth, input expression stratum, and FDR threshold. Each user-defined scenario is repeated for 100 times, and these metrics are computed and averaged to generate empirical estimations.

Given a MeRIP-seq dataset in .bam files, various experimental scenarios (such as sample size, sequencing depth, FDR threshold, etc.), and a chosen differential methylation testing method, magpie generates evaluation results for each proposed study design. Functions incorporated in magpie allow users to export these results in an .xlsx file, and to visualize them through line plots. Users have the option to provide small pilot data, which could include only several chromosomes. We would estimate major parameters from these pilot data, to guide larger-scale simulations for power evaluation for future experimental designs. Alternatively, when pilot MeRIP-seq datasets are unavailable or unattainable, the function quickPower can produce power evaluation results within seconds. This is achieved by directly extracting our in-house evaluation results based on three 8 re .0 w oP .60

0 C .1 D F .5 0

B D public N6-methyladenosine datasets on GEO as the pilot data (Niu et al., 2013; Schwartz et al., 2014; Barbieri et al., 2017) . Our package also comes with a vignette that provides thorough instructions and examples of its applications in differential analysis experimental design on N6-methyladenosine.

Under simulation settings outlined in Supplementary Materials S1, we next examine the relationship between sample size and power in DMR detection, given that determining sample size is a primary objective in our method. We adopt sample sizes of 2, 3, 5, 7, and 10 per group, and nominal FDR values of 0.05, 0.1, 0.15, and 0.2, both of which are common choices in MeRIP-seq experiments. The empirical results for major metrics are shown in Figure 2. Grouped by sample size and nominal FDR level, power, targeted power, FDC, and FDR averaged over 100 simulations are shown in Figure 2AD. For a fixed sample size, metrics like power, FDC, and FDR diminish under lower FDR thresholds. This occurs as lower FDR values lead to greater stringency, which in turn reduces false positives. The power will drop, as expected, when using stringent FDRs. As sample size increases, these differences become smaller, particularly for statistical power (Figure 2A). Here, power remains consistently high across all FDR levels with 7 and 10 replicates per group. This highlights the benefit of using a larger sample size that helps detect DMRs with limited effect sizes, where a type II error would often occur when the sample size is small. Such trend is observed consistently when using different pilot data (Figure S1). At the same time, these results give researchers the knowledge to optimize the sample size based on their budgets. Using Figure 2A as an exam

Mean Input: (30, 51] FDR Threshold = 0.05 8 .

0 8 .

0 C D F .4 0 0 . 0 ple, a power around 0.8 is achieved with 7 samples per group, and sample size of 7 is considered large in current MeRIP-seq studies. The benefit of expanding the sample size to 10 is marginal, but the associated costs could be significantly higher.

It is useful for researchers to understand the effects of heterogeneity in baseline expression levels in DMR detection. In MeRIP-seq data, the basal expression level is represented by input control read counts, thus we stratify power metrics by input control ranges. Six strata are obtained based on following quantiles of mean input counts: stratum 1 (0%-10%), stratum 2 (10%-30%), stratum 3 (30%-50%), stratum 4 (50%-70%), stratum 5 (70%-90%), and stratum 6 (90%-100%). At a nominal FDR of 0.05, the average targeted power and FDC for the six strata are shown in Figure 3A, B. Overall, reduced targeted power is observed in the lower strata, a trend that is more evident when sample sizes are small. This is expected, as true differences in low-expressed regions are often obscured by noise, making DMRs harder to detect. A limited sample size further exacerbates this issue. This suggests the potential benefits of increasing the sequencing depth, particularly when biological replicates are limited and more samples are hard to obtain. Here, relatively low strata will enjoy the benefit of more drastic power improvement. Interestingly, higher FDCs are reported in the upper strata, suggesting that more false positives are detected per true discovery in these highly expressed regions. However, this trend diminishes with increasing sample size. Given that these metrics were computed across various simulation scenarios, we further explore the variability of the results, using visualizations within a specific stratum and sample size in Figure 3C-F. With elevated sample sizes, there is reduced variability in both targeted power and FDC (Figure 3C, D). This is not surprising since it is more likely to capture the true dispersion with more replicates, leading to more consistent power estimates. However, this trend is not observed across the strata at a fixed sample size, suggesting the benefit of increasing sample size over sequencing depth for more reliable inferences. Under a fixed sample size (Figure 3E, F), an upward trend is observed across the strata for both targeted power and FDC. This trend aligns with the observations in Figure 3A, B, though some variability is evident. A heatmap panel is also available to illustrate the stratified results (Figure S3).

It is worth noting that the targeted power and FDC presented in the Results section are computed for DMRs with odds ratios (OR) exceeding Δ = 2, using TRESS. To evaluate the fluctuations of these two metrics across various effect sizes (OR), sample sizes and DMR detection methods, we also consider Δ values of 1.5, 2, 4, 6, 8, and 10, for both TRESS and exomePeak2. In Figure 4, the targeted power and FDC are plotted against the sample size and are grouped by odds ratio thresholds. At all sample sizes, there is an increased targeted power (Figure 4A, C) and a higher FDC to identify DMRs

C D A B 0 .

1 rew .80 oP .60 ted .4 e 0 ragT .02 0 . 0 2 0 1 8 C 6 FD 4 2

2 0

0 2 3 5 7 3 5 7

10 10 2

Sample Size with a larger odds ratio. Specifically, for FDC (Figure 4B, D), substantially higher values are observed among DMRs with exceptionally large odds ratios (Δ = 8, 10). This indicates that detecting DMRs with these large ORs might lead to a significant increase in false positives. These patterns hold true for both TRESS and exomePeak2. While both methods show improvements in targeted power with added replicates across all odds ratio thresholds, a discrepancy is noted for FDC that it tends to increase with larger sample sizes when using exomePeak2. This discrepancy, however, is not universally observed when applying our proposed framework to different pilot data sets (Figure S2). These findings highlight the importance of utilizing the users’ designated DMR detection methods during power calculation to ensure accurate estimations.

As shown previously, sequencing depth is another critical factor in MeRIP-seq study design. Building upon our analysis of sequencing coverage strata, here we examine another aspect of sequencing depth by introducing a “depth factor". This is a relative ratio to reflects the effect to enlarge or down-sample the sequencing coverage of the pilot data. As illustrated in Figure 5A, the targeted power rises with increased sequencing depth in all sample sizes. The incremental gain from increasing sequencing depth diminishes at high depths or large sample sizes, but benefits the small sample size the most. In Figure 5B for FDC, a similar pattern is observed as in the stratified analysis: FDCs increase with sequencing depth, but stabilize in scenarios with larger sample sizes. We also provide an integrated visualization in Figure 5C, presenting targeted power and FDC in the same panel, aiding users in understanding the tradeoffs between them. Researchers could consult similar figures, generated by magpie using their own pi5