A Hybrid HMM Approach for the Dynamics of DNA Methylation

TL;DR

This paper introduces a hybrid hidden Markov model to analyze the rapid and complex dynamics of DNA methylation and demethylation in cells, integrating deterministic and stochastic processes.

Contribution

The paper presents a novel hybrid HMM that models both deterministic and random events in DNA methylation dynamics, tested on mouse embryonic stem cell data.

Findings

Predicted methylation changes over short time scales.

Estimated efficiencies of methylation and demethylation steps.

Validated model on time-resolved data from stem cells.

Abstract

The understanding of mechanisms that control epigenetic changes is an important research area in modern functional biology. Epigenetic modifications such as DNA methylation are in general very stable over many cell divisions. DNA methylation can however be subject to specific and fast changes over a short time scale even in non-dividing (i.e. not-replicating) cells. Such dynamic DNA methylation changes are caused by a combination of active demethylation and de novo methylation processes which have not been investigated in integrated models. Here we present a hybrid (hidden) Markov model to describe the cycle of methylation and demethylation over (short) time scales. Our hybrid model decribes several molecular events either happening at deterministic points (i.e. describing mechanisms that occur only during cell division) and other events occurring at random time points. We test our model on mouse embryonic stem cells using time-resolved data. We predict methylation changes and estimate the efficiencies of the different modification steps related to DNA methylation and demethylation.

Tables

Table 1. Table 1: Cell division matrix 𝐃 𝐃 \mathbf{D} .

	$uu$	$um$	$uh$	$uf$	$mu$	$mm$	$mh$	$mf$	$hu$	$hm$	$hh$	$hf$	$fu$	$fm$	$fh$	$ff$
$uu$	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$um$	1/2	1/2	0	0	0	0	0	0	0	0	0	0	0	0	0	0
$uh$	1/2	0	1/2	0	0	0	0	0	0	0	0	0	0	0	0	0
$uf$	1/2	0	0	1/2	0	0	0	0	0	0	0	0	0	0	0	0
$mu$	1/2	0	0	0	1/2	0	0	0	0	0	0	0	0	0	0	0
$mm$	0	1/2	0	0	1/2	0	0	0	0	0	0	0	0	0	0	0
$mh$	0	0	1/2	0	1/2	0	0	0	0	0	0	0	0	0	0	0
$mf$	0	0	0	1/2	1/2	0	0	0	0	0	0	0	0	0	0	0
$hu$	1/2	0	0	0	0	0	0	0	1/2	0	0	0	0	0	0	0
$hm$	0	1/2	0	0	0	0	0	0	1/2	0	0	0	0	0	0	0
$hh$	0	0	1/2	0	0	0	0	0	1/2	0	0	0	0	0	0	0
$hf$	0	0	0	1/2	0	0	0	0	1/2	0	0	0	0	0	0	0
$fu$	1/2	0	0	0	0	0	0	0	0	0	0	0	1/2	0	0	0
$fm$	0	1/2	0	0	0	0	0	0	0	0	0	0	1/2	0	0	0
$fh$	0	0	1/2	0	0	0	0	0	0	0	0	0	1/2	0	0	0
$ff$	0	0	0	1/2	0	0	0	0	0	0	0	0	1/2	0	0	0

Table 2. Table 2: Conversion errors for CpGs, where the rates for single cytosines are defined in Fig.  4 . We define x ¯ := 1 − x assign ¯ 𝑥 1 𝑥 \bar{x}:=1-x . Note that for MAB-Seq we also define j := μ ​ d + ( 1 − μ ) ​ ( 1 − c ) assign 𝑗 𝜇 𝑑 1 𝜇 1 𝑐 j:=\mu d+(1-\mu)(1-c) .

	bisulfite seq. (BS)				ox. bisulfite seq. (oxBS)				MAB-Seq
	TT	TC	CT	CC	TT	TC	CT	CC	TT	TC	CT	CC
$uu$	$c^2$	$c\cdot \overline{c}$	$c\cdot \overline{c}$	${\overline{c}}^2$	$c^2$	$c\cdot \overline{c}$	$c\cdot \overline{c}$	${\overline{c}}^2$	${\overline{j}}^2$	$j\cdot \overline{j}$	$j\cdot \overline{j}$	$j^2$
$um$	$c\cdot \overline{d}$	$c\cdot d$	$\overline{c}\cdot \overline{d}$	$\overline{c}\cdot d$	$c\cdot \overline{d}$	$c\cdot d$	$\overline{c}\cdot \overline{d}$	$\overline{c}\cdot d$	$\overline{j}\cdot \overline{d}$	$\overline{j}\cdot d$	$j\cdot \overline{d}$	$j\cdot d$
$uh$	$c\cdot \overline{e}$	$c\cdot e$	$\overline{c}\cdot \overline{e}$	$\overline{c}\cdot e$	$c\cdot f$	$c\cdot \overline{f}$	$\overline{c}\cdot f$	$\overline{c}\cdot \overline{f}$	$\overline{j}\cdot \overline{e}$	$\overline{j}\cdot e$	$j\cdot \overline{e}$	$j\cdot e$
$uf$	$c\cdot g$	$c\cdot \overline{g}$	$\overline{c}\cdot g$	$\overline{c}\cdot \overline{g}$	$c\cdot g$	$c\cdot \overline{g}$	$\overline{c}\cdot g$	$\overline{c}\cdot \overline{g}$	$\overline{j}\cdot g$	$\overline{j}\cdot \overline{g}$	$j\cdot g$	$j\cdot \overline{g}$
$mu$	$c\cdot \overline{d}$	$\overline{c}\cdot \overline{d}$	$c\cdot d$	$\overline{c}\cdot d$	$c\cdot \overline{d}$	$\overline{c}\cdot \overline{d}$	$c\cdot d$	$\overline{c}\cdot d$	$\overline{j}\cdot \overline{d}$	$j\cdot \overline{d}$	$\overline{j}\cdot d$	$j\cdot d$
$mm$	${\overline{d}}^2$	$d\cdot \overline{d}$	$d\cdot \overline{d}$	$d^2$	${\overline{d}}^2$	$d\cdot \overline{d}$	$d\cdot \overline{d}$	$d^2$	${\overline{d}}^2$	$d\cdot \overline{d}$	$d\cdot \overline{d}$	$d^2$
$mh$	$\overline{d}\cdot \overline{e}$	$\overline{d}\cdot e$	$d\cdot \overline{e}$	$d\cdot e$	$\overline{d}\cdot f$	$\overline{d}\cdot \overline{f}$	$d\cdot f$	$d\cdot \overline{f}$	$\overline{d}\cdot \overline{e}$	$\overline{d}\cdot e$	$d\cdot \overline{e}$	$d\cdot e$
$mf$	$\overline{d}\cdot g$	$\overline{d}\cdot \overline{g}$	$d\cdot g$	$d\cdot \overline{g}$	$\overline{d}\cdot g$	$\overline{d}\cdot \overline{g}$	$d\cdot g$	$d\cdot \overline{g}$	$\overline{d}\cdot g$	$\overline{d}\cdot \overline{g}$	$d\cdot g$	$d\cdot \overline{g}$
$hu$	$c\cdot \overline{e}$	$\overline{c}\cdot \overline{e}$	$c\cdot e$	$\overline{c}\cdot e$	$c\cdot f$	$\overline{c}\cdot f$	$c\cdot \overline{f}$	$\overline{c}\cdot \overline{f}$	$\overline{j}\cdot \overline{e}$	$j\cdot \overline{e}$	$\overline{j}\cdot e$	$j\cdot e$
$hm$	$\overline{d}\cdot \overline{e}$	$d\cdot \overline{e}$	$\overline{d}\cdot e$	$d\cdot e$	$\overline{d}\cdot f$	$d\cdot f$	$\overline{d}\cdot \overline{f}$	$d\cdot \overline{f}$	$\overline{d}\cdot \overline{e}$	$d\cdot \overline{e}$	$\overline{d}\cdot e$	$d\cdot e$
$hh$	${\overline{e}}^2$	$e\cdot \overline{e}$	$e\cdot \overline{e}$	$e^2$	$f^2$	$f\cdot \overline{f}$	$f\cdot \overline{f}$	${\overline{f}}^2$	${\overline{e}}^2$	$e\cdot \overline{e}$	$e\cdot \overline{e}$	$e^2$
$hf$	$\overline{e}\cdot g$	$\overline{e}\cdot \overline{g}$	$e\cdot g$	$e\cdot \overline{g}$	$f\cdot g$	$f\cdot \overline{g}$	$\overline{f}\cdot g$	$\overline{f}\cdot \overline{g}$	$\overline{e}\cdot g$	$\overline{e}\cdot \overline{g}$	$e\cdot g$	$e\cdot \overline{g}$
$fu$	$c\cdot g$	$\overline{c}\cdot g$	$c\cdot \overline{g}$	$\overline{c}\cdot \overline{g}$	$c\cdot g$	$\overline{c}\cdot g$	$c\cdot \overline{g}$	$\overline{c}\cdot \overline{g}$	$\overline{j}\cdot g$	$j\cdot g$	$\overline{j}\cdot \overline{g}$	$j\cdot \overline{g}$
$fm$	$\overline{d}\cdot g$	$d\cdot g$	$\overline{d}\cdot \overline{g}$	$d\cdot \overline{g}$	$\overline{d}\cdot g$	$d\cdot g$	$\overline{d}\cdot \overline{g}$	$d\cdot \overline{g}$	$\overline{d}\cdot g$	$d\cdot g$	$\overline{d}\cdot \overline{g}$	$d\cdot \overline{g}$
$fh$	$\overline{e}\cdot g$	$e\cdot g$	$\overline{e}\cdot \overline{g}$	$e\cdot \overline{g}$	$f\cdot g$	$\overline{f}\cdot g$	$f\cdot \overline{g}$	$\overline{f}\cdot \overline{g}$	$\overline{e}\cdot g$	$e\cdot g$	$\overline{e}\cdot \overline{g}$	$e\cdot \overline{g}$
$ff$	$g^2$	$g\cdot \overline{g}$	$g\cdot \overline{g}$	${\overline{g}}^2$	$g^2$	$g\cdot \overline{g}$	$g\cdot \overline{g}$	${\overline{g}}^2$	$g^2$	$g\cdot \overline{g}$	$g\cdot \overline{g}$	${\overline{g}}^2$

Equations

$$\mathbf{P}(t)=\mathbf{D}\cdot \mathbf{M}(t).$$

$$\pi (t_i+\Delta t)=\pi (t_i)\cdot \mathbf{P}(t_i),$$

$$\frac{d}{dt}\pi (t)=\pi (t)\cdot \mathbf{Q}(t).$$

$$\underset{t\to t_i^{-}}{\lim }\pi (t)=\pi (t_i^{-})\ne \pi (t_i^{-})\cdot \mathbf{P}(t_i^{-})=\pi (t_i^{+})=\underset{t\to t_i^{+}}{\lim }\pi (t).$$

$$r(t):={\alpha }_r+{\beta }_r\cdot t,$$

$${\pi }_e(t)=\pi (t)\cdot E_e,$$

$$\pi (0{)}^{\ast }={\text{argmax}}_{\pi (0)}{\mathcal{L}}_1(\pi (0)),$$

$${\mathcal{L}}_1(\pi (0))=\prod _{e\in E}\prod _{j\in {\mathcal{S}}_{obs}}{\pi }_e(j,0{)}^{n_e(j,0)}.$$

$${\vec{v}}^{\ast }={\text{argmax}}_{\vec{v}}{\mathcal{L}}_2(\vec{v}),$$

$${\mathcal{L}}_2(\vec{v})=\prod _{e\in E}\prod _{t\in T_{obs}\backslash \{0\}}\prod _{j\in {\mathcal{S}}_{obs}}{\pi }_e(j,t{)}^{n_e(j,t)}.$$

Full text

Here we present a hybrid (hidden) Markov model to describe the cycle of methylation and demethylation over (short) time scales. Our hybrid model decribes several molecular events either happening at deterministic points (i.e. describing mechanisms that occur only during cell division) and other events occurring at random time points. We test our model on mouse embryonic stem cells using time-resolved data. We predict methylation changes and estimate the efficiencies of the different modification steps related to DNA methylation and demethylation.

1 Introduction

All cells of a multi-cellular organism share the same DNA sequence, yet, depending on location and cell type, display distinct cellular programs as a result of controlled gene expression. Hence, the expression of genes is regulated by epigenetic factors such as DNA methylation. In mammals, the methylation of DNA is restricted to the C5 position of cytosine (C) and mostly appears in a CpG di-nucleotide sequence [6, 7]. The palindromic nature of CpG positions provides a symmetry which, after DNA replication, allows the stable inheritance of the methylation “signal”. Methylation of C to 5-methyl cytosine (5mC) is catalysed by a certain enzyme family, the DNA methyltransferases (Dnmts) [2, 28]. Three conserved and catalytic active family members are associated with the methylation of DNA. Dnmt1 is mainly responsible for maintenance methylation after DNA replication, i.e. the enzyme mainly reestablishes the methylation pattern on the newly synthesized daughter strand [15] according to the inherited information of the parental DNA strand. The enzymes Dnmt3a and Dnmt3b perform de novo methylation, where new methyl groups are added on unmethylated Cs [27]. However, there is evidence that this separation of tasks is not definite and that all Dnmts may carry out all tasks to a certain extent [1].

Once established, 5mC can be further modified by oxidation to 5-hydroxymethyl cytosine (5hmC), which can again be oxidized to 5-formyl cytosine (5fC). Then, 5fC can eventually be converted to 5-carboxy cytosine (5caC). All these processes are carried out by the ten-eleven translocation (Tet) enzymes [16, 29]. A considerable level of 5hmC can be found in many cells types and its occurrence has been connected to gene regulation as well as genome wide loss of DNA methylation [12, 19]. In contrast, 5fC and 5caC are far less abundant and their particular functions remain more illusive. Nevertheless, studies suggest that both oxidized cytosine variants function as intermediates during enzymatic removal of 5mC from the DNA [12].

In general, DNA methylation can be removed in two ways. First, after DNA replication, the absence or blocking of maintenance methylation will result in a passive DNA methylation loss with each cell division (passive demethylation). Second, generated 5fC or 5caC is enzymatically removed from the DNA and subsequently replaced by unmodified cytosine (active demethylation) [4, 13, 14, 24].

While DNA replication and the associated maintenance methylation happens only once per cell division cycle, de novo methylation and the modification processes of 5mC via oxidation, as well as active demethylation, may happen at arbitrary time points. With purely discrete hidden Markov models (HMM), as used in [1] and [9], it is difficult to describe multiple instances of de novo methylation or other modification events for a given CpG during one cell division cycle. Here, we introduce a hybrid HMM to describe the dynamics of active demethylation. It distinguishes between events at fixed time points (cell division and maintenance methylation) and events at random time points (de novo methylation, oxidations, active demethylation). By applying our model to a data set of a single copy gene from mouse embryonic stem cells, we were able to accurately predict the frequency of the observable CpG states and the levels of the hidden states, which correspond to the different modified forms of C. Furthermore, we show how to estimate the enzymatic reaction efficiencies using a maximum likelihood approach.

Compared to previous models for describing DNA methylation, we here propose an approach that takes into account both, the rather fixed and deterministic timing of cell division and the random nature of the enzymatic processes. In this way, we are able to present a model that realistically describes the dynamics of DNA methylation and improves our understanding of active demethylation.

The paper is organized as follows: In Section 2 we give the necessary biological and mathematical background, i.e. we explain passive and active demethylation in more detail, describe the model and explain the parameter estimation procedure. In Section 3 the results are discussed and in Section 4 we conclude our findings.

2 Model

2.1 Passive and active demethylation

One can distinguish between two different ways of losing methylation at cytosines: After cell division a new (daughter) strand of the DNA is synthesized. Initially, all cytosines of the daughter strand are unmethylated, while the methylation states of cytosines on the parental strand remain unchanged. Maintenance methylation, which happens during the replication process, is used to re-establish the methylation pattern at the newly synthesized strand. However, the absence or inhibition of maintenance methylation causes a loss of 5mC with each replication step. This DNA replication dependent loss has been termed passive demethylation.

For active demethylation it is assumed that oxidation of 5mC to 5fC or 5caC via 5hmC and a subsequent enzymatic removal of the oxidative cytosine from the DNA occurs (cf. Fig 1). Since de novo methylation as well as active demethylation are replication independent, the loop depicted in Fig. 1 can be traversed multiple times within one cell cycle. In measurements, we see Cs in all stages, i.e., either methylated, oxidized (5hmC/5fC/5caC) or unmodified.

2.2 Hybrid Markov model

In this section, we present a model that describes the state changes of a single CpG over time. It can be seen as a hybrid extension of previous discrete-time models [1, 9]. The (hidden) states of a CpG correspond to the set of all pairs of the four possibilities in Fig. 1, i.e., $\{u,m,h,f\}^{2}$, because it contains a C on both strands of the DNA. We split the transitions of our model into transitions that occur at fixed times and those that occur at random times. This results in a mixture of a discrete time Markov chain (DTMC) and a continuous time Markov chain (CTMC). In the following we will refer to the events or transitions that occur at fixed time points as discrete part of the model, while we refer to the other events or transitions at the random time points as continuous part of the model.

We assume that cells divide after a fixed time interval (usually every 24 hours). Hence, these events correspond to deterministic transitions at fixed times. During cell division one strand is kept as it is (parental strand) and all methylation states and its modifications remain unchanged, while one DNA strand is newly synthesized (daughter strand) and therefore contains only unmethylated cytosine. Consequently, after cell divison a CpG that was modified on both sides becomes a CpG that is unmodified on one side. Since the parental strand is chosen at random, the probability for each of the two successor states (corresponding to the state of the CpG in the two daughter cells) is 0.5. The full transition probability matrix $\mathbf{D}$ for cell division is shown in Tab. 1. Note that the cell division matrix is time independent.

Maintenance methylation, i.e. methylation events that occur on hemimethylated CpGs to reestablish methylation patterns, is known to be linked to the replication fork [22]. We therefore consider maintenance to occur together with the cell division at the same fixed time points. Hence, cell division and maintenance can be described by a (discrete-time) Markov chain whose transition probability matrix $\mathbf{P}(t)$ is defined in the sequel in Eq. (1). Maintenance may happen at the daughter strand if there is a methylated C on the parental strand, i.e. on hemimethylated CpGs ($um$ or $mu$), with probability $\mu_{m}(t)$. Since it is reasonable to assume that hemihydroxylated CpGs ($uh$ or $hu$) have different properties compared to hemimethylated CpGs in terms of maintaining existing methylation patterns, we describe the maintenance probability for hemihydroxylated CpGs as follows. Let $p$ be the probability that 5hmC is recognized as unmethylated by maintenance enzymes, i.e., the enzyme will not perform maintenance of a hemihydroxylated CpG. Then, the maintenance probability of such a CpG is given by $\bar{p}\mu_{m}(t)$, where $\bar{p}=1-p$. Note that there is no equivalent probability for $uf$ or $fu$, i.e. we assume that CpGs with 5fC or 5caC at one strand are not maintained [17]. The transition probability matrix for maintenance events is illustrated in Fig. 2, where we omitted the time dependency of $\mu_{m}(t)$. Whenever no transition is possible, i.e. there is only a self loop for this state (omitted in Fig. 2) the corresponding diagonal entry in the matrix is 1. Where a transition is possible we set the corresponding (off-diagonal) entry in the matrix to its transition probability and the diagonal entry to 1 minus its transition probability. All other entries in the matrix are 0.

One discrete step of the corresponding DTMC corresponds to one cell division, including maintenance methylation. Hence, its transition probability matrix is defined as

[TABLE]

Every other event may occur an arbitrary (unknown) number of times between two cell divisions at random time points and will be described by a continuous-time Markov jump process. These events are de novo methylation ($u\rightarrow m$) with rate $\mu_{d}(t)$, hydroxylation ($m\rightarrow h$) with rate $\eta(t)$, formylation ($h\rightarrow f$) with rate $\phi(t)$ and active demethylation ($f\rightarrow u$) with rate $\delta(t)$ (cf. Fig 1). Note that all these events may happen on both strands, independent of the state on the complementary strand. All possible reactions are shown in Fig. 3. From these transitions the infinitesimal generator matrix $\mathbf{Q}(t)$ of the jump process can easily be inferred. For the off-diagonal elements, we set the entries to the respective reaction rate if a reaction is possible between two states, indicated by a colored arrow in Fig. 3, and to [math] if no reaction is possible. The diagonal elements are then given by the negative sum of the off-diagonal elements of the respective row.

In order to describe the time evolution let us first define the set of time points $T_{d}=\{t_{1},t_{2},\ldots,t_{n}\}$, at which cell division and maintenance occur. We assume that there are in total $n$ of these events. At these time points $t_{i}$ the time evolution of the probability distribution of the states is given by

[TABLE]

where $\Delta t$ is the duration of cell division and maintenance. After cell division and maintenance the other events take place at random time points until the time point for the next cell division is reached. During this interval $[t_{i}+\Delta t,t_{i+1}]$ the time evolution for $\pi(t)$ is obtained by solving the differential equation

[TABLE]

Note that since cell division and maintenance methylation occur at a time interval that is much shorter than the time between two cell divisions, we assume here that these events occur instantaneously, i.e. we let $\Delta t\to 0$, which leads to a jump of the distribution at $t_{i}$. It holds that the left-hand and right-hand limit do not coincide (except for the special case of only $uu$ at time $t_{i}$), i.e.

[TABLE]

To resolve the ambiguity we set the value of $\pi$ at time $t_{i}$ to $\pi(t_{i}):=\lim_{t\rightarrow t_{i}^{+}}\pi(t)$, which means that we assume that at time $t_{i}$ the cell division and maintenance methylation have already happened. Intuitively, we can obtain the solution for $\pi(t)$ numerically by alternating between multiplication with $\mathbf{P}$ (Eq. (2)) and integration of Eq. (3).

2.3 Efficiencies

As already discussed in [9] the methylation or modification rates may change over time, i.e. with the assumption of constant rates the behavior of the system can not be captured correctly. We therefore introduce time dependent rates, which we call efficiencies. The simplest way of making the efficiencies time dependent is to impose a linear form. Let $r\in\{\mu_{m},\mu_{d},\eta,\phi,\delta\}$ be one of the reaction rates. For each reaction we then define

[TABLE]

where $\alpha_{r}$ and $\beta_{r}$ are some new parameters to characterize the linear efficiency function. To further compress the notation we write these new parameters into a vector $\vec{v}_{r}=(\alpha_{r},\beta_{r})$.

Note that in order to ensure identifiability of the parameters we have to introduce certain constrains. Obviously since $\mu_{m}(t)$ and $p$ are probabilities, they have to be bound between zero and one, i.e. $0\leq\mu_{m}(t)\leq 1$ for all $t\in[0,t_{max}]$ and $0\leq p\leq 1$. The other efficiencies have to be bound by some upper limit $ub$ since otherwise the methylation cycle may become arbitrarily fast and the optimization algorithm runs into identifiability problems. We therefore require $0\leq r(t)\leq ub$ for all $t\in[0,t_{max}]$, where $r(t)\neq\mu_{m}(t)$. Since typical average turnover times $\mathbb{E}[T_{turnover}]=\mu_{d}^{-1}+\eta^{-1}+\phi^{-1}+\delta^{-1}$ of the methylation cycle are in the order of $75$ to $120$ minutes in certain promotors of human cells [18, 25], a viable choice would be, for example, $ub=12$, i.e. each modification occurs on average not more frequent than 12 times per hour (not faster than every $5$ minutes).

2.4 Conversion errors

The actual state of a CpG can not be directly observed. We therefore have to estimate the hidden states from sequencing experiments. Since the different modifications of C might lead to the same observable states we perform three different kinds of sequencing experiments: All sequencing strategies share a bisulfite treatment, which usually converts C and its oxidized variants 5fC and 5caC, summarized to 5fC∗, to uracil. Additionally, in bisulfite sequencing (BS) both 5mC and 5hmC remain unconverted while in oxidative bisulfite sequencing (oxBS) only 5mC is retained as C [3]. A combination of both methods can therefore be used to estimate the amount of 5hmC. In M. SssI assisted bisulfite sequencing (MAB-Seq) at first all Cs in a CpG context are methylated and afterwards BS is applied. Thus, if all conversions would happen without any errors, only 5fC∗ would be converted to T. In order to capture the methylation pattern of CpG position at both complementary DNA strands, the distinct chemical treatments were combined with hairpin sequencing [10, 11, 21]. Regular reactions with their respective probabilities are marked with solid black arrows in Fig. 4, while the possible false reactions are depicted with dashed red arrows. Since every CpG consists of two Cs (one on each strand) with independent conversion errors, we get the conversion error for CpGs by multiplying the individual conversion errors. A complete overview of all possible combinations for each of the three methods is shown in Tab. 2.

2.5 Parameter estimation

Recall that the set of hidden states is $\mathcal{S}=\{uu,um,uh,uf,mu,mm,mh,mf,hu,$ $hm,hh,hf,fu,fm,fh,ff\}$. We now define a hidden Markov model (HMM) based on the model presented in Section 2.2. As set of observable states we define $\mathcal{S}_{obs}=\{TT,TC,$ $CT,CC\}$, i.e. we use the results of the sequencing experiments (cf. Fig. 4) on both strands. The conversion errors define the corresponding emission probabilities. We also define $n_{e}(j,t)$ as the number of times that state $j\in S_{obs}$ has been observed during independent measurements of sequencing method $e\in E:=\{\text{BS, oxBS, MAB-Seq}\}$. The probability distribution over all observable states for experiment $e$ is denoted by $\pi_{e}(t)$, the probability of a state $j\in S_{obs}$ by $\pi_{e}(j,t)$ and in a similar fashion we denote $\pi(i,t)$ with $i\in S$ for the hidden states, with probability distribution $\pi(t)$. The observable and hidden states for all times $t$ are connected via

[TABLE]

where $E_{e}$ is the emission matrix for sequencing method $e$ and is listed in Tab. 2 for each of the three methods.

Our goal is to estimate the efficiencies for the different methylation events given our hybrid HMM and data from the three different experiments at different time points $t\in T_{obs}$ via a maximum likelihood estimator (MLE). Since an initial distribution over the hidden states, which can not directly be observed, is needed in order to initialize the model, we have to employ the MLE twice: First we estimate the initial distribution $\pi(0)$ over the hidden states by maximizing

[TABLE]

under the constrain that $\sum_{i\in S}\pi(i,0)=1$. The likelihood $\mathcal{L}_{1}(\pi(0))$ is defined as

[TABLE]

Note that Eq. (8) is independent of the parameters. Given an initial distribution over the hidden states we can now run our model and apply the MLE

[TABLE]

a second time in order to estimate the efficiencies, where

[TABLE]

The vector $\vec{v}=(\vec{v}_{\mu_{m}},\vec{v}_{\mu_{d}},\vec{v}_{\eta},\vec{v}_{\phi},\vec{v}_{\delta},p)$ contains all unknown parameters for all efficiencies and the probability $p$ of considering a hydroxylated cytosine as unmethylated. Note that applying the MLE twice and independently leads only to an approximation of the true most likely explanation, since the estimated initial distribution may not lead to the same result in the parameter estimation, as if it would all be done in one estimation. However, we choose this approach to reduce the computational complexity of the optimization. In order to estimate the standard deviations of the estimated parameters we use the observed Fisher information matrix [5]. The Fisher information is defined as $\mathcal{J}(\vec{v}^{*})=-\mathcal{H}(\vec{v}^{*})$, where $\vec{v}^{*}$ is the maximum likelihood estimate and $\mathcal{H}(\vec{v})=\nabla\nabla^{T}\log\mathcal{L}_{2}(\vec{v})$ the Hessian matrix of the log-likelihood. The expected Fisher information is then given by $\mathcal{I}(\vec{v}^{*})=\mathbb{E}[\mathcal{J}(\vec{v}^{*})]$ and its inverse forms a lower bound for the covariance matrix. Thus, we can approximate the standard deviation of all estimated parameters by $\sigma(\vec{v}^{*})=\sqrt{\text{Var}(\vec{v}^{*})}\approx\sqrt{\text{diag}(-\mathcal{H}^{-1}(\vec{v}^{*}))}$.

The implementation of the hybrid HMM and its analysis as explained above has been integrated into the latest beta version of the H(O)TA tool [20]. H(O)TA provides results for individual CpGs and also an aggregated profile across all analyzed CpGs.

3 Results

In the following, we will discuss the results after applying our model to data derived from a short region at the single copy gene Afp (alpha fetoprotein), which contains 5 CpGs. More precisely, we followed the DNA methylation changes during the adjustment of mouse embryonic stem cells (mESCs) towards 2i medium after previous long time cultivation under Serum/LIF conditions [8, 9]. The Serum/ LIF-to-2i shift is a common model system which induces genome wide demethylation in mESCs including the Afp locus.

Here, the individual and aggregated H(O)TA results show a very similar behavior. Hence, we only present the aggregated results in Fig. 5.

Fig. 5 (a) - (c) shows a comparison of the actual measurements (solid lines) and the results from the model (dashed lines) of the time evolution for the four observable states for all three sequencing models. The predictions and measurements are in good agreement. In BS and oxBS over time the frequency of TT increases, while the frequency of CC decreases. Moreover, the frequencies of TC and CT show a temporal increase. For MAB-Seq the frequencies of all observable states remain quite constant, except for some small changes in early times. This behavior can be explained by the changes over time in the efficiencies, which are shown in Fig. 5 (d). The maintenance probability $\mu_{m}$ remains constant at about $0.8$, as well as the probability $p$ of considering 5hmC as unmethylated, which is constant by definition. The estimation of $p$ gives $p=1$, i.e. 5hmC is always considered to be unmethylated and will not be recognized by Dnmt1 after replication which means that 5hmC leads to an impairment of Dnmt1 activity and a passive loss of DNA methylation with each replication. The de novo efficiency $\mu_{d}$ decreases over time, while the hydroxylation and formylation efficiency $\eta$ and $\phi$ increase. The demethylation efficiency $\delta$ is always 0. Note, however, that standard deviation becomes very large for later time points due to insufficient data. The efficiencies explain the behavior of the frequencies of the observable states and is even more evident for the hidden states shown in Fig. 5 (e): Over time the probability of being fully or hemimethylated decreases, while the probability of being unmethylated, hydroxylated or formylated increases. A more detailed look into hydroxylated and formylated states is shown in Fig. 5 (f) and (g). Note that the combination of hydroxylation and formylation is only shown in one of the two subplots, namely in (g). The observed increase in 5hmC and 5fC/5caC is in accordance with the high oxidation efficiency of Tets in form of hydroxylation and formylation. Previously, we showed using a purely discrete HMM, that the presence of 5hmC leads to a block of Dnmt1 activity after replication, which we also observe in the present hybrid model [9]. However, we now also observe an increase in higher oxidized cytosine variants, namely 5fC and 5caC which equally prevents methylation by Dnmt1. Thus, we reason that the impact of Tet mediated oxidation of 5mC on DNA demethylation in the investigated system plays a much more important role than previously suggested [26].

Considering the rather rapid periodic events of de novo methylation and active demethylation, the chosen measurement time points are not ideal. There is only information available at the end of each cell division cycle (one division within 24h), i.e., no information is given for the times between two cell divisions. With measurements at time points between two cell divisions we would be able to distinguish if a CpG is in a certain state because it initially was, or because it ran through the full cycle, possibly multiple times. We would like to emphasize that with better data, i.e. with sufficiently many measurements between cell divisions, it will be possible to estimate all reaction efficiencies with better confidence.

4 Conclusion

We proposed a hybrid hidden Markov model which is able to successfully describe both, events such as cell division and maintenance methylation that occur at fixed times, and events that occur at random times, such as de novo methylation, oxidizations and active demethylation, according to a continuous-time Markov jump process. To the best of our knowledge, this is the first model that describes the dynamics of active demethylation, i.e., the active removal of the methyl group through several enzymatic steps. We applied our model to data from mouse embryonic stem cells, which undergo a gradually loss of DNA methylation over time. We were able to accurately predict the frequency of the observable states and the levels of the hidden states in all cases. We were also able to predict the enzymatic reaction efficiencies based on a linear assumption for their time behavior.

As future work we plan to apply our model to more informative data such that all efficiencies of the active demethylation cycle can be estimated with better confidence. Moreover, we plan to allow different functional forms for the efficiencies since the linear form is not flexible enough and does not allow to capture complex behaviors. A suitable choice could be splines of different degrees and with a different number of knots. However, in this case it is also necessary to perform model selection in order to prevent overfitting. Another possible extension could be to investigate potential neighborhood dependencies of the modified cytosines [23].