N ⁶ -methyl-adenosine (m⁶A) is the most abundant modification among 100 types of identified RNA modifications in eukaryotic mRNA/lncRNA [1, 2]. Even though m⁶A was found existing in mammalian mRNAs in as early as 1970s [3], its biological relevance remains unclear due to the difficulties in identifying global m⁶A sites in mRNA [4]. In 2013, the m⁶A demethylase Fat mass and obesity associated protein (FTO) was first discovered [5], to be able to reverse the m⁶A modification in mRNA and it revived our interests of studying m⁶A in mRNA. To date, ALKBH5 is identified as another demethylase [6] and the methyltransferase like 3/14 (METTL3/METTL14) and Wilms’ tumor 1-assoicating protein (WTAP) are discovered to be subunits of the m⁶A methyltransferase complex [7, 8]. All these findings provide strong evidences to show that m⁶A is a dynamic modification and suggest that it may play a critical role in exerting post-transcriptional functions in mRNA metabolism [9–11].

These new wave of breakthroughs cannot be achieved without the recent development of MeRIP-seq [12, 13], which was successfully developed to reveal the transcriptome-wide distribution of m⁶A in human and mouse cells. In this essay, mRNA is first chemically fragmented into approximately 100-nucleotide (nt) long before immunoprecipitation with anti-m⁶A antibody. Then, the immunoprecipitated (IPed) methylated mRNA fragments and the un-immunoprecipited input control mRNA fragments are subjected to high-throughput sequencing [14]. The sequenced IP and input reads are aligned to the transcriptome and reads enrichment of IP out of the combined reads in IP and input samples are examined to predict to loci of methylation sites and infer the degree of methylation. We have previously developed exomePeak [15, 16] and HEPeak [17], two algorithms for detecting m⁶A peaks in MeRIP-seq. Although MeRIP-seq and subsequent computational peak-calling analysis provide an accurate landscape of m⁶A methylation in transcriptome, the complete mechanisms of this methylation still remains unclear. Just like gene expression where co-expression might suggest co-regulation or similar gene functions, sites with similar methylation degree could be related to similar methylation mechanisms. Therefore, there is a need to develop algorithms to uncover co-methylation pattern in MeRIP-seq data. In this paper, we model the methylation degrees of m⁶A peaks as a mixture of the Beta-binomial distributions and propose an expectation-maximization based clustering algorithm to uncover the co-methylation patterns.

In this section, we first describe the proposed generative model to define m⁶A peak clusters and then derive the Expectation-Maximization algorithm for the inference. In the end, we discuss a Bayesian Information Criterion (BIC) [18] for selecting the optimal number of m⁶A peak clusters.

The proposed graphical model for clustering RNA methylation peaks

The proposed graphical model for clustering of m⁶A peaks in MeRIP-seq data is shown in Fig. 1. Suppose we have identified a set of N m⁶A peaks, by using peak-calling software such as exomePeak or HEPeak. The goal is to cluster these peaks according to their methylation degree, which is defined as IP reads count divided by the total count of IP and control reads. For the n _th m⁶A peak, let Z _n  ∈ {1, 2,.., K} denote the index of the particular methylation cluster that n-th peak belongs to, with K representing the total number of clusters, then Z _n follows a discrete distribution

$$ P\left({Z}_n\Big|\boldsymbol{\uppi} \right)={\displaystyle \prod_{k=1}^K{\pi_k}^{I\left({Z}_n=k\right)}} $$

(1)

where π _k is the unknown probability that an m⁶A peak belongs to cluster k, where ∑ _K π _k  = 1 and I(⋅) is the indicator function. Also, let the observed reads count in the n _th peak of the m _th IP replicate sample be X _m,n and that of the m _th input replicate denote as Y _mn . Under the assumption that reads count follows a Poisson distribution, the reads count X _mn given the total reads account T _mn  = X _mn  + Y _mn can be shown to follow a Binomial distribution

$$ P\left({X_m}_n\Big|{p}_n,{Z}_n\right)={{\displaystyle \prod_{k=1}^K\left(\left(\begin{array}{c}\hfill {T}_{mn}\hfill \\ {}\hfill {X}_{mn}\hfill \end{array}\right){p_n}^{X_{mn}}{\left(1-{p}_n\right)}^{Y_{mn}}\right)}}^{I\left({Z}_n=k\right)} $$

(2)

where is p _n represents unknown methylation degree at the n _th Peak of the m _th replicate. In order to model the variance of the replicates for the n _th peak, given cluster assignment Z _n , p _n is assumed to follow the Beta distribution

$$ P\left(p\Big|{Z}_n\right)={\displaystyle \prod_{k=1}^K Beta{\left({\alpha}_k,{\beta}_k\right)}^{I\left({Z}_n=k\right)}} $$

(3)

Therefore, after integrating the variable \( {p}_n \), \( {X}_{mn} \) follows a mixture of Beta-binomial distribution

$$ \begin{array}{l}P\left({X}_{mn}\Big|{Z}_n;\boldsymbol{\upalpha}, \boldsymbol{\upbeta} \right)={\displaystyle \sum_{p_n}P\left({X_m}_n\Big|{p}_n,{Z}_n\right)P\left({p}_n\Big|\boldsymbol{\upalpha}, \boldsymbol{\upbeta} \right)}\\ {}={{\displaystyle \prod_{k=1}^K\left(C\bullet \frac{\varGamma \left({X}_{mn}+{\alpha}_k\right)\varGamma \left({Y}_{mn}+{\beta}_k\right)}{\varGamma \left({T}_{mn}+{\alpha}_k+{\beta}_k\right)}\frac{\varGamma \left({\alpha}_k+{\beta}_k\right)}{\varGamma \left({\alpha}_k\right)\varGamma \left({\beta}_k\right)}\right)}}^{I\left\{{Z}_n=k\right\}}\end{array} $$

(4)

where α = [α ₁, α ₂,.., α _K ] ^T , β = [β ₁, β ₂,.., β _K ] ^T are the unknown parameters of Beta distribution and C is the normalization constant. Thus, by considering the N m⁶A peaks in M replicates, the joint distribution is

$$ P\left(\mathbf{X},\mathbf{Z}\Big|\boldsymbol{\upalpha}, \boldsymbol{\upbeta}, \boldsymbol{\uppi} \right)={\displaystyle \prod_{n=1}^N{\displaystyle \prod_{m=1}^M{\displaystyle \prod_{k=1}^K{\left({\pi}_kBB\left({X}_{mn}\Big|{Z}_n\right)\right)}^{I\left\{{Z}_n=k\right\}}}}} $$

(5)

where BB(X _mn |Z _n ) represents formula (3). Then, the log-likelihood of the observed data can be expressed as

$$ \begin{array}{l}l= \lg P\left(\mathbf{X}\Big|\boldsymbol{\upalpha}, \boldsymbol{\upbeta}, \boldsymbol{\uppi} \right)= \lg {\displaystyle \sum_{\mathbf{Z}}P\left(\mathbf{X},\mathbf{Z}\Big|\boldsymbol{\upalpha}, \boldsymbol{\upbeta}, \boldsymbol{\uppi} \right)}\\ {}={\displaystyle \sum_{n=1}^N{\displaystyle \sum_{m=1}^M \lg {\displaystyle \sum_{k=1}^K{\pi}_kBB\left({X}_{mn}\Big|{Z}_n;\boldsymbol{\upalpha}, \boldsymbol{\upbeta} \right)}}}\end{array} $$

(6)

where Z = [Z ₁, Z ₂, …, Z _N ] ^T , X = [X _1, ^T X _2,…, ^T X _N ^T ] ^T and \( {\boldsymbol{X}}_{\boldsymbol{n}}={\left[{X}_{1n},{X}_{2n},\dots, {X}_{Mn}\right]}^T \). The goal of inference is to predict the cluster index \( {Z}_n \) for all the peaks and estimate the unknown model parameters \( \boldsymbol{\theta} =\left[\boldsymbol{\alpha}, \boldsymbol{\beta}, \boldsymbol{\pi} \right] \). Next, we first discuss the maximum likelihood solution for parameter inference, based which an EM algorithm is introduced afterwards to perform model parameters inference and cluster assignment jointly.

In this paper, a novel graphical model based methylation peak clustering algorithm, was developed for discovering the patterns in methylation degrees of m⁶A peaks in the MeRIP-seq data. The peak cluster is modelled as the mixture Beta-binomial distribution, where the Beta distribution can model the variance of the methylation degree across sample replicates. The evaluation on both simulation and real MeRIP-seq datasets demonstrates the accuracy and robustness of our model. In addition, our algorithm successfully uncovered a unique and novel pattern for m⁶A peak cluster, providing a new lead for understanding the mechanisms and functions of m⁶A methylation.

BIC, Bayesian Information Criterion; CDS, Coding DNA sequence; EM, Expectation of maximum likelihood method; FDR, False discovery rate; MeRIP-seq, Methylated RNA Immunoprecipatation combined with RNA sequencing; UTR, Untranslated region