Dependency Modeling

Modeling site-site dependencies with heterogeneous HMM and bidirectional LSTM

Code · Thesis

This work uses a heterogeneous HMM and a bidirectional LSTM to model how methylation states co-vary across nearby CpG sites on individual sequencing reads. The trained models can serve as conditional generative models that produce realistic read-level methylation patterns while preserving both site-level averages and local dependencies.

Motivation

Site-Level Summaries Lose Read-Level Structure

Bisulfite sequencing records binary methylation states on individual DNA fragments. Conventional analysis compresses these reads into one methylation level per CpG, calculated as the fraction of methylated reads covering that site. This aggregation discards information about which methylation states co-occur on the same read.

Two distinct read-level methylation patterns with identical methylation levels at every site
Two read-level patterns with the same site-level methylation profile.

Every site in the example has a methylation level of 0.5, yet Pattern 1 is mixed while Pattern 2 contains coordinated blocks. The summary retains each marginal probability $P(Y_t=1)$ but discards joint information such as $P(Y_t,Y_{t+1})$, so matching the site-level averages does not guarantee realistic read-level patterns.

Why Independent Sampling Is Not Enough

A common baseline models each CpG as an independent Bernoulli variable parameterized by its site-level methylation probability. This preserves the marginals by construction but removes co-methylation—the tendency of nearby CpGs on the same read to share a state.

Box plots comparing dependency statistics in real WGBS data and independent Bernoulli samples
Dependency statistics for real WGBS reads and independent Bernoulli samples.

The modeling objective therefore has two parts:

  • the average generated state at each site should match its observed methylation level
  • the generated joint patterns should reproduce how dependency changes with genomic distance and sequence context

Models

Diagrams of independent Bernoulli, heterogeneous HMM, and bidirectional LSTM methylation models
Independent Bernoulli, distance-aware HMM, and feature-rich BiLSTM models.

Heterogeneous Hidden Markov Model

The model separates the true methylation state $z_t$ from the observed state $y_t$. Its transition matrix changes with the genomic distance $d_t$ between adjacent CpG sites:

\[A(d_t)=A_1+\phi(d_t)A_2, \qquad \phi(d_t)=\operatorname{sigmoid}(w_0+w_1d_t),\quad w_1<0.\]

$A_1$ describes the distance-independent transition pattern, while $A_2$ supplies a distance-dependent adjustment. Nearby sites are more likely to retain the same state; as distance grows, $\phi(d_t)$ decreases and the transition matrix approaches $A_1$.

The emission model maps the hidden biological state to the state observed after bisulfite conversion and sequencing:

\[B= \begin{bmatrix} 1-p_3 & p_3 \\ p_4 & 1-p_4 \end{bmatrix}.\]

Here, $p_3$ and $p_4$ summarize incomplete or inappropriate bisulfite conversion together with sequencing error. This separation keeps biological dependency in the transition model and experimental noise in the emission model.

Bidirectional LSTM

The model uses a feature vector $x_t$ containing the methylation level, genomic distance, local sequence context, and flanking bases at each site. It processes the sequence in both directions:

\[\overrightarrow{h}_t =\operatorname{LSTM}_{\mathrm{forward}}(x_t,\overrightarrow{h}_{t-1}), \qquad \overleftarrow{h}_t =\operatorname{LSTM}_{\mathrm{backward}}(x_t,\overleftarrow{h}_{t+1}).\]

The two representations are concatenated and mapped to a methylation probability:

\[h_t=[\overrightarrow{h}_t;\overleftarrow{h}_t], \qquad \ell_t=W_{\mathrm{fc}}h_t+b_{\mathrm{fc}}, \qquad o_t=\sigma(\ell_t).\]

This allows the model to learn bidirectional patterns beyond the HMM’s first-order, distance-only assumption.

Stochastic Pattern Generation

At generation time, temperature scaling and Gumbel noise introduce controlled variation into the output probabilities:

\[g_t=-\log(-\log u_t),\quad u_t\sim\operatorname{Uniform}(0,1), \qquad \widetilde{p}_t=\sigma\!\left(\frac{\ell_t}{\tau}+\gamma g_t\right).\]

The temperature $\tau$ controls randomness, while $\gamma$ determines the contribution of the noise. Repeated draws yield diverse read-level methylation patterns from the same input features while retaining the dependencies learned by the BiLSTM.

Implementation

Heterogeneous HMM Estimation

The parameters $\theta=(p_1,p_2,w_0,w_1,p_3,p_4)$ are estimated by maximum likelihood. Treating the $N$ reads as independent, the observed-data likelihood is

\[\mathcal{L}(\theta) =\prod_{r=1}^{N}P(Y_r\mid D_r,\theta), \qquad P(Y_r\mid D_r,\theta) =\sum_{Z_r}P(Y_r,Z_r\mid D_r,\theta).\]

For a single read, the complete-data likelihood factorizes into initial-state, emission, and distance-dependent transition terms:

\[P(Y_r,Z_r\mid D_r,\theta) =P(z_{r,1})P(y_{r,1}\mid z_{r,1}) \prod_{t=2}^{T_r} P(z_{r,t}\mid z_{r,t-1},d_{r,t})P(y_{r,t}\mid z_{r,t}).\]

Because the observed-data likelihood sums over all possible hidden patterns $Z_r$, the model uses an Expectation-Maximization algorithm to maximize the expected complete-data log-likelihood:

\[Q(\theta\mid\theta^{(m)}) =\mathbb{E}_{Z\mid Y,D,\theta^{(m)}} \left[\log P(Y,Z\mid D,\theta)\right].\]

The E-step uses a forward-backward pass to compute hidden-state and transition posteriors. The M-step updates emission probabilities in closed form and optimizes the distance-dependent transition parameters with L-BFGS-B. Numerical clipping and multiple initializations improve stability; the highest-likelihood fit is retained.

BiLSTM Training Objective

The BiLSTM was implemented in PyTorch with padded mini-batches and masks for variable-length reads. Its custom objective is

\[L_{\text{total}}=\alpha L_{\text{BCE}}+\beta L_{\text{MSE}}, \qquad \alpha=\beta=1.\]

$L_{\text{BCE}}$ learns the read-level binary states, while $L_{\text{MSE}}$ aligns the average of stochastic predictions with the target site-level methylation scores. Training uses an 80/20 split, dropout, gradient clipping, and held-out loss monitoring over 50 epochs.

Evaluation

A controlled synthetic-data study evaluates two aspects of the learned models: whether their generated sequences preserve the target marginal methylation profile and whether they reproduce the known distance-dependent inter-site relationship.

Comparison of true and predicted marginal methylation distributions across four models on synthetic data
Marginal preservation across models on synthetic data.
Probability of adjacent sites sharing the same state across synthetic data and four model predictions
Distance-dependent relationships in synthetic data and model predictions.

On synthetic data, the BiLSTM reproduced the dependency-distance correlation ($r=-0.92$ versus $r=-0.95$ in the ground truth) while maintaining strong marginal alignment ($r=0.91$, MAE $=0.137$).

Summary

Together, the plots show the central trade-off: matching site-level averages is not the same as preserving their joint structure. The heterogeneous HMM offers an interpretable model of distance-dependent relationships, while the BiLSTM provides the flexibility to incorporate additional features and capture more complex patterns. Beyond describing these dependencies, the fitted models can generate realistic read-level methylation patterns for downstream method development and evaluation.