**4. Discussion**

Results from our simulation study indicate that BFM is a valid approach to detect DMRs when considering ordinal group responses, since the calculated Bayes factors were very large for simulated DMRs, and close to 1 for non-DMRs. The real data analysis based on the CLL data also demonstrated that BFM is a valid method that is able to detect DMRs with methylation rates increasing (or decreasing) as disease severity increases.

In addition to being able to account for ordering of group responses, BFM also has an advantage of allowing for heterogeneity of methylation effects across CpG sites by modeling the methylation rates with a prior. Methods such as the SSM pools information across variants in a region, assuming that each CpG sites in the region have the same methylation rates.

BFM with mixed-effect regression, not only can allow for covariates, but also the correlation between CpG sites. It takes advantage of the flexibility of the Bayesian framework, including the use of prior information when available as well as computational convenience, and uses distributions such as the multivariate normal to incorporate the correlation structure with inverse Wishart distribution as the prior for the correlation matrix.

One disadvantage of the BFM is that it assumes that methylation rates of CpG sites within each moving window are independent of those outside of the window, while the SSM accounts for the correlation along the whole genome.

BFM used a moving window to help decide the location and length of DMRs. But practically, it is very difficult to know the exact length of DMRs. This limitation is very common in statistical genetics, not only for detecting DMRs, but also for detecting rare variants [63]. Cross validation or bootstrap approaches might help determine the window sizes. It could be possible to develop other methods, for example, using genes and promoters instead of moving windows, along with the BFM to detect DMRs.

As described by George and Laud [64], the Bayes factor used in the context of testing hypotheses is a meaningful measure of evidence because it is a reasonably approximate factor by which the odds are increased by the data. With default priors that are essentially flat over a wide range of the relevant parameter space, the approach is similar to the likelihood-based inference. However, direct comparison between methods such as BFM based on the Bayes factor with frequentist approaches should be done with caution, as the Bayes factor classification for decision process is not a precise calibration, but rather a descriptive statement about the standards of evidence. Our proposed method is rather exploratory in nature, leading to a ranked list of sites for follow up for formal confirmation.

We developed the BFM, focusing only on DNA methylation data. However, large-scale cancer genomics projects such as TCGA (The Cancer Genome Atlas Research Network) are currently generating multiple layers of genomics data for early tumor, including DNA copy number, methylation, and mRNA expression. Similar statistical methods for integrated analysis and systematic modeling of these genomics data deserve further attention.

**Author Contributions:** Conceptualization, F.D., H.X. and V.G.; methodology, F.D., H.X, D.R., S.G. and V.G.; formal analysis, F.D.; writing—original draft preparation, F.D.; writing—H.X. and V.G.; data curation, H.S.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors extend our special thanks and appreciation to one of the reviewers who provided a very extensive review of the manuscript with several insightful and constructive comments, which helped us improve the manuscript substantially.

**Conflicts of Interest:** The authors declare no conflict of interest.
