### HEPeak pipeline

To address the aforementioned MeRIP-seq issues, HEPeak includes several high-throughput sequencing tools in its pipeline. First, HEPeak utilizes TopHat [18] to align fragmented mRNA reads to the reference transcriptome, allowing short reads to span exon-exon junctions. Next, SAM-tools [19] is applied to exclude the multi-mapping reads and index alignment results. After these pre-processing steps, HEPeak performs HMM-based peak calling on the exons of each gene, where the introns are excluded, to identify the genomic locus of methylation sites. The output result of HEPeak is in BED format, which can be visualized together with input alignments in IGV2.1 [20].

### Exome-based peak calling

The goal of peak calling in MeRIP-seq is to detect regions in transcripts where the read counts in the IP sample is more "enriched" than those in the control sample. Just as with ExomePeak, our previously developed peak calling algorithm for MeRIP-seq, HEPeak performs the peak calling on connected exons of a specific gene, a clear contrast to genome-based ChIP-seq peak calling methods, such as MACS [21]. This projection of genome onto transcriptome effectively circumvents the difficulty due to the ambiguity of isoforms' assignment but it still preserves the convenience of gene-based annotation, making biological interpretation of the prediction straightforward.

### The definition of HMM for MeRIP-seq data

Given a particular mRNA (RefSeq gene), its concatenated exons are first divided into *N* mutually connected bins, whose size is selected as the read length *L*. With respect to the *n*_{
th
} bin, the unknown hidden methylation status is denoted as *z*_{
n
} ∈ {1, 2} where 1 represents unmethylation and 2 otherwise. Since a peak likely spans multiple bins, we assume that the methylation status *z*_{
n
} follows a first order Markov chain, whose transition matrix *A* contains entries defined as

{A}_{jk}=P\left({z}_{n}=k|{z}_{n-1}=j\right),\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}\phantom{\rule{2.77695pt}{0ex}}j,k\in \left\{1,2\right\}

(1)

where *A*_{
jk
} denotes the probability for the latent variable switching from the status *j* at the (*n* - 1)_{
th
} bin to the status *k* at the *n*_{
th
} bin. Here *j*, *k* is the indicator of the hidden state. Additionally, we assume that the initial probability *P*(*z*_{l} = 1) = *π* and *P*(*z*_{l} = 2) = 1 - *π*.

Next, let *x*_{
n
} denote the read counts in the IP sample and *y*_{
n
} the counts in the control sample, both for bin *n*. We assume that, given the methylation status *z*_{
n
}, these read counts follow the Poisson distribution defined as

P\left({x}_{n}|{t}_{n}\right)=Pois\left({M}_{IP}{\lambda}_{IP,{z}_{n}}\right)

(2)

P\left({y}_{n}|{t}_{n}\right)=Pois\left({M}_{ctrl}{\lambda}_{ctrl}\right)

(3)

where *M*_{
IP
} and *M*_{
ctrl
} are the total reads (sequencing depth) in the IP and the control samples, respectively and {\lambda}_{IP,{z}_{n}} for *z*_{
n
} = 1, or 2 and *λ*_{
ctrl
} are the normalized Poisson rates, respectively. It is worthwhile pointing out that {\lambda}_{IP,{z}_{n}} switches according to the status of *z*_{
n
}; on the contrary, *λ*_{
ctrl
} stays the same.

It would be intuitive next to define the relationship between the Poisson rates for the methylated and unmethylated in the IP and the control sample, respectively. However, unlike in ChIP-seq, where this relationship is mostly defined only for the IP sample, defining the relationship for both the IP and the control is non-trivial and model complexity also needs to be assessed to avoid potential difficulties in subsequent inference. To this end, we transform the formulation by observing that, given (2) and (3), the conditional probability of observing *x*_{
n
} in the IP given the total reads in the control as *t*_{
n
} = *x*_{
n
} + *y*_{
n
} follows the binomial distribution

P\left({x}_{n}|{z}_{n},{t}_{n}\right)=Bino\left({t}_{n},{p}_{{z}_{n}}\right)

(4)

where

{p}_{{z}_{n}}=\frac{{M}_{IP}{\lambda}_{IP,{z}_{n}}}{{M}_{ctrl}{\lambda}_{ctrl}+{M}_{IP}{\lambda}_{IP,{z}_{n}}}.

(5)

Note that {p}_{{z}_{n}} for *z*_{
n
} = 1 (or 2) can be considered as the percentage of the mean IP read counts in the combined read counts of the IP and control samples for a bin, when it is unmethylated (or methylated). The distribution (4) effectively combines the reads in the IP and control samples under one model. As such, instead of using (2) and (3), we define (4) as the emission probability of the proposed HMM and work with {p}_{{z}_{n}} directly. Doing so avoids modelling and inferring the potentially complex relationships between the rates. Given **X** = {*x*_{1}, *x*_{2}, *x*_{3},..., *x*_{
N
}}, a set of reads for *N* bins and **Z** = {*z*_{1}, *z*_{2}, *z*_{3},..., *z*_{
N
}}, the sequence of methylation, we use *γ*(*z*_{n,k}) to denote the marginal posterior distribution of a latent variable *z*_{
n
}at state *k*, and *ε*(*z*_{n-1}, *z*_{
n
}) to denote the joint posterior distribution of two successive latent variables, so that

\gamma \left({z}_{n,k}\right)=p\left({z}_{n}=k|X,\theta \right)

(6)

\epsilon \left({{z}_{n-1}}_{,j},{z}_{n,k}\right)=p\left({z}_{n-1}=j,{z}_{n}=k|X,\theta \right).

(7)

Here, the parameter is defined as \theta =\left\{{A}_{k,j}\forall k\forall j;\pi ;{p}_{k}\forall k\right\}. Then, the log likelihood for the proposed HMM chain can be expressed as

\begin{array}{c}Q={E}_{z}\left[\mathsf{\text{ln}}P\left(X,Z|\theta \right)\right]={\displaystyle \sum _{k=1}^{2}}\gamma \left({z}_{1,k}\right)\mathsf{\text{ln}}{\pi}_{k}\hfill \\ \phantom{\rule{2.77695pt}{0ex}}+{\displaystyle \sum _{n=1}^{N}}{\displaystyle \sum _{j=1}^{2}}\gamma \left({z}_{n,k}\right)\mathsf{\text{ln}}P\left({x}_{n}|{z}_{n,k}\right)+{\displaystyle \sum _{n=2}^{N}}{\displaystyle \sum _{j=1}^{2}}{\displaystyle \sum _{k=1}^{2}}\epsilon \left({z}_{n-1,j}{z}_{n,k}\right)\mathsf{\text{ln}}{A}_{jk}\hfill \end{array}

(8)

We call this new formulation HEPeak or Hidden Markov Model (HMM)-based Exome Peak finding. The graphical model of HEPeak formulation is shown in Figure 2A. Compared with ExomePeak, HEPeak considers the correlation of the reads between adjacent bins and more accurately models the behaviour of methylated reads in MeRIP-Seq (Figure 2B).

### The EM solution

Given HEPeak, the goal is to call peaks, i.e., predict *z*_{
n
}∀*n*, and at the same time estimate the model parameters: *θ*. To this end, we developed an Expected-Maximization (EM) solution, which performs peak calling and parameter estimation in an iterative fashion. We provide the steps of the EM algorithm in the following. The detailed derivation is included in appendix.

At the *m*_{
th
} iteration, proceed as follows.

**E step**: Given parameter *θ*^{(m-1)}, estimated at the *m*-1 step, calculate the posterior distribution of the latent variable *P*(**Z**|**X**, *θ*^{(m-1)}).

\gamma \left({z}_{n,k}\right)=p\left({z}_{n}=k|X,{\theta}^{\left(m-1\right)}\right)

(9)

**M step**: Compute and update *π*^{(m)}, *A*_{
jk
}^{(m) }and *p*_{
k
}^{(m) }for all *j*, *k* as

\pi =\frac{\gamma \left({z}_{11}\right)}{{\displaystyle \sum _{j=0}^{1}}\gamma \left({z}_{1j}\right)}

(10)

{A}_{jk}=\frac{{\displaystyle \sum _{n=2}^{N}}\epsilon \left({z}_{n-1,j},{z}_{n,k}\right)}{{\displaystyle \sum _{l=0}^{1}}{\displaystyle \sum _{n=2}^{N}}\epsilon \left({z}_{n-1,j},{z}_{n,l}\right)}

(11)

{p}_{k}=\frac{{\displaystyle \sum _{n=1}^{N}}\gamma \left({z}_{n,k}\right){X}_{n}}{{\displaystyle \sum _{n=1}^{N}}\gamma \left({z}_{nk}\right)\left({X}_{n}+{Y}_{n}\right)}

(12)

After the EM iteration converges, the model parameter *θ* can be obtained. Given the estimated *θ*, the Viterbi algorithm is applied to maximize the joint likelihood in (8) to obtain the maximum *a posteriori* (MAP) estimate of the methylation status *z*_{
n
}.

### Peak region detection

In order to evaluate the statistical significance of the putative peak regions predicted by the Viterbi algorithm, the log odds ratio of the posterior for the peak state (*z*_{
n
} = 2) over the posterior for the background state (*z*_{
n
} = 1) can be computed as follows

\mathsf{\text{PeakScore}}\left({z}_{n}\right)=\mathsf{\text{log}}\frac{p\left({z}_{n}=2|X\right)}{p\left({z}_{n}=1|X\right)}

(13)

Briefly, this log-transformed scoring method [22–24] tries to utilize the posterior probability of each bin to assess the confidence of the potential peak region. The potential peak region is defined as consecutive bins predicted by the Viterbi and its PeakScore is calculated as the averaged PeakScores for all the combined bins. Next, PeakScore is assumed to follow a Gaussian distribution with mean (*mean*(PeakScore)) and standard deviation(*std*(PeakScore)) [24], estimated from all the bins. Then, after performing the z transform of PeakScores, a one-sided test for significance of the potential peak region can be conducted and p-value can calculated. Then, the Benjamini-Hochberg method [25] is utilized to correct the multiple testing and compute the False Discovery Rate (FDR).