Beyond similarity assessment: Selecting the optimal model for sequence alignment via the Factorized Asymptotic Bayesian algorithm

TL;DR

This paper introduces a new model selection method for Pair Hidden Markov Models in sequence alignment, improving accuracy by choosing optimal hidden state configurations using the Factorized Information Criteria.

Contribution

The paper presents a novel model selection approach for PHMMs that effectively determines the number of hidden states, enhancing alignment accuracy in bioinformatics.

Findings

The method accurately selects models with higher posterior probability.

Simulation results show improved alignment accuracy.

Application to DNA datasets demonstrates the selection of more complex models.

Abstract

Pair Hidden Markov Models (PHMMs) are probabilistic models used for pairwise sequence alignment, a quintessential problem in bioinformatics. PHMMs include three types of hidden states: match, insertion and deletion. Most previous studies have used one or two hidden states for each PHMM state type. However, few studies have examined the number of states suitable for representing sequence data or improving alignment accuracy.We developed a novel method to select superior models (including the number of hidden states) for PHMM. Our method selects models with the highest posterior probability using Factorized Information Criteria (FIC), which is widely utilised in model selection for probabilistic models with hidden variables. Our simulations indicated this method has excellent model selection capabilities with slightly improved alignment accuracy. We applied our method to DNA datasets from 5 and 28 species, ultimately selecting more complex models than those used in previous studies.

Tables

Table 1. Table 1: Model sizes of hand-crafted models. The 1st, 2nd and 3rd values in each triplet are the numbers of match states ( K M subscript 𝐾 𝑀 K_{M} ), X 𝑋 X -insertion states ( K X subscript 𝐾 𝑋 K_{X} ) and Y 𝑌 Y -insertion states ( K Y subscript 𝐾 𝑌 K_{Y} ), respectively. Model names including the term ”imb” indicate an imbalanced model with unequal values of K X subscript 𝐾 𝑋 K_{X} and K Y subscript 𝐾 𝑌 K_{Y} . See Supplementary Figures S1–S7 for details of the parameters.

small	med	large	imb	imb_large	huge	imb_huge
(1, 1, 1)	(1, 2, 2)	(1, 4, 4)	(1, 2, 1)	(1, 4, 2)	(1, 6, 6)	(1, 6, 3)

Table 2. Table 2: Fraction of precise model selections from multiple estimations for data produced with different model sizes. We ran the model selection algorithm 10 times for each combination of sample number and model size (shown in Table 1 ) with random initial parameters. Each table entry indicates the fraction of runs in which the correct optimal model was selected. The bold-face font indicates precise prediction of the model size, that is when the method selected the model with the largest FIC correctly.

# samples	precise model selection for different model sizes
	small	med	large	imb	imb_large	huge	imb_huge
100	10/10	0/10	0/10	10/10	0/10	0/10	0/10
200	10/10	0/10	0/10	9/10	0/10	0/10	0/10
300	10/10	0/10	0/10	4/10	7/10	0/10	0/10
400	10/10	0/10	0/10	9/10	10/10	0/10	10/10
500	10/10	0/10	0/10	7/10	9/10	0/10	5/10
600	10/10	5/10	0/10	9/10	10/10	5/10	9/10
700	10/10	10/10	3/10	5/10	10/10	10/10	10/10
800	10/10	10/10	10/10	7/10	10/10	10/10	10/10
900	9/10	10/10	10/10	8/10	10/10	10/10	10/10
1000	9/10	8/10	10/10	7/10	10/10	10/10	10/10

Table 3. Table S1: Insertion f1 scores

	Trained models
Simulated model	small	med	large	imb	imb_large	huge	imb_huge	fab
small	0.5725	0.5764	0.5676	0.5689	0.5715	0.5653	0.5677	0.5716
med	0.7414	0.7401	0.7375	0.7446	0.7371	0.7337	0.7355	0.7404
large	0.6365	0.6631	0.6685	0.6506	0.6653	0.6640	0.6695	0.6701
imb	0.6456	0.6775	0.6733	0.6767	0.6710	0.6731	0.6733	0.6778
imb_large	0.6105	0.6422	0.6428	0.6330	0.6427	0.6390	0.6393	0.6412
huge	0.8623	0.8736	0.8792	0.8684	0.8759	0.8825	0.8792	0.8821
imb_huge	0.9000	0.9025	0.9111	0.9027	0.9104	0.9089	0.9094	0.9091
average	0.7380	0.7494	0.7502	0.7462	0.7500	0.7480	0.7495	0.7514

Equations

$$\displaystyle\ln p(\bm{X},\bm{Z}|\bm{\Pi})=\sum_{n=1}^{N}\bigg{[}\ln p(\bm{z}^{n}_{\text{in}}|\bm{\alpha})+\sum_{t=0}^{T_{X}^{n}}\sum_{u=0}^{T_{Y}^{n}}\Big{(}\ln p(\bm{z}^{n}_{tu}|pa(\bm{z}^{n}_{tu}),\bm{\beta})+\ln p(x^{n}_{t},y^{n}_{u}|\bm{z}^{n}_{tu},\bm{\phi})\Big{)}\bigg{]}$$

$$p(x_t,y_u|z_{tuk}^n=1)={\psi }_{tuk}^n=\{\begin{array}{l}{\varphi }_k(x_t,y_u)\text{if}\mathcal{S}(k)=M\\ {\varphi }_k(x_t,-)\text{ if}\mathcal{S}(k)=X\\ {\varphi }_k(-,y_u)\text{ if}\mathcal{S}(k)=Y\end{array}$$

$$\ln p(\mathbf{X},\mathbf{Z}|\mathbf{\Pi })$$

$$\displaystyle=\sum_{n=1}^{N}\sum_{k=1}^{K}\bigg{[}z^{n}_{d_{xk},d_{yk},k}\ln\alpha_{k}+\sum_{t=0}^{T_{X}^{n}}\sum_{u=0}^{T_{Y}^{n}}\Big{(}\sum_{l=1}^{K}z^{n}_{(t-d_{xk}),(u-d_{yk}),k}z^{n}_{tul}\ln\beta_{kl}+z^{n}_{tuk}\ln\psi_{tuk}^{n}\Big{)}\bigg{]}$$

$$(d_{xk},d_{yk})=\{\begin{array}{l}(1,1)\text{if}\mathcal{S}(k)=M\\ (1,0)\text{if}\mathcal{S}(k)=X\\ (0,1)\text{if}\mathcal{S}(k)=Y.\end{array}$$

$${\mathcal{M}}^{\ast }=\underset{\mathcal{M}}{\mathrm{argmax}}\ln p(\mathbf{X}|\mathcal{M})$$

$$\displaystyle FIC(\mathcal{M})=\mathbb{E}_{q^{**}}\big{[}\ln p(\bm{X},\bm{Z}|\bar{\bm{\Pi}},\mathcal{M})-\frac{1}{2}\ln|F_{\bar{\bm{\Xi}}}|\big{]}-\frac{1}{2}D_{\bm{\Pi}}\ln N+\mathcal{H}(q^{**}).$$

$$F_{{\mathbf{\beta }}_k}=\frac{1}{N}\sum _{n=1}^N\sum _{t=0}^{T_X}\sum _{u=0}^{T_Y}pa({\mathbf{z}}_{tu}^n{)}_k{\nabla }_{\mathbf{\beta }}^2\ln p_k({\mathbf{z}}_{tu}^n|{\mathbf{\beta }}_k)$$

$$F_{{\mathbf{\varphi }}_k}=\frac{1}{N}\sum _{n=1}^N\sum _{t=0}^{T_X}\sum _{u=0}^{T_Y}{z}_{tuk}^n{\nabla }_{\mathbf{\varphi }}^2\ln p(x_t^n,y_u^n|{\mathbf{z}}_{tu}^n,{\mathbf{\varphi }}_k),$$

$$\ln |F_{{\mathbf{\beta }}_k}|=D_{{\mathbf{\beta }}_k}\ln \frac{{\zeta }_k^{\text{trans}}(\mathbf{Z})}{N}+\mathcal{O}(1)$$

$$\ln |F_{{\mathbf{\varphi }}_k}|=D_{{\mathbf{\varphi }}_k}\ln \frac{{\zeta }_k^{\text{emit}}(\mathbf{Z})}{N}+\mathcal{O}(1)$$

$$\{\begin{array}{ll}{\zeta }_k^{\text{trans}}(\mathbf{Z})& ={\sum }_{n=1}^N{\sum }_{t=0}^{T_X}{\sum }_{u=0}^{T_Y}pa({\mathbf{z}}_{tu}^n{)}_k\\ & ={\sum }_{n=1}^N{\sum }_{t=0,u=0}^{T_x^n,T_y^n}{z}_{tuk}^n-{\sum }_{n=1}^N{z}_{T_x^n,T_y^n}^n.\\ {\zeta }_k^{\text{emit}}(\mathbf{Z})& ={\sum }_{n=1}^N{\sum }_{t=0,u=0}^{T_x^n,T_y^n}{z}_{tuk}^n.\end{array}$$

$$\displaystyle FIC_{\textit{PHMM}}(\mathcal{M})=\mathbb{E}_{q^{**}}\bigg{[}\ln p(\bm{X},\bm{Z}|\bar{\bm{\Pi}})-\sum_{k=1}^{K}\frac{D_{\beta_{k}}}{2}\ln\frac{\zeta^{\text{trans}}_{k}(\bm{Z})}{N}-\sum_{k=1}^{K}\frac{D_{\phi_{k}}}{2}\ln\frac{\zeta^{\text{emit}}_{k}(\bm{Z})}{N}\bigg{]}-\frac{D_{\Pi}}{2}\log N+\mathcal{H}(q^{**}).$$

$$FI{C}_{\text{PHMM}}(\mathcal{M})$$

$$\displaystyle=\mathbb{E}_{q}\bigg{[}\log p(\bm{X},\bm{Z}|\bm{\Pi})+\sum_{n=1}^{N}\sum_{t=0,u=0}^{T^{n}_{x},T^{n}_{y}}z^{n}_{tuk}\log\delta_{tuk}\bigg{]}-\frac{D_{\alpha}}{2}\log N$$

$$\displaystyle\qquad-\sum_{k=1}^{K}\frac{D_{\beta_{k}}}{2}\log\big{(}\sum_{n,t,u}\zeta_{k}^{\text{trans}}(\tilde{\bm{Z}})-1\big{)}-\sum_{k=1}^{K}\frac{D_{\phi_{k}}}{2}\log\big{(}\zeta_{k}^{\text{emit}}(\tilde{\bm{Z}})-1\big{)}+\mathcal{H}(q)\quad$$

$$\displaystyle\text{where}\quad\delta_{tuk}=\left\{\begin{aligned} &\exp\bigg{(}-\frac{D_{\phi_{k}}}{2\zeta_{k}^{\text{emit}}(\tilde{\bm{Z}})}\bigg{)}\enskip\text{if}\enskip t=T_{X}\enskip\text{and}\enskip u=T_{Y}\\ &\exp\bigg{(}-\frac{D_{\phi_{k}}}{2\zeta_{k}^{\text{emit}}(\tilde{\bm{Z}})}-\frac{D_{\beta_{k}}}{2\zeta_{k}^{\text{trans}}(\tilde{\bm{Z}})}\bigg{)}\enskip\text{otherwise}.\end{aligned}\right.$$

$$f_{tuk}^n$$

$$b_{tuk}^n$$

$$q^{\ast }(z_{tuk}^n)$$

$$q^{\ast }(z_{t-d_{xk},u-d_{xy},j}^n,z_{tuk}^n)$$

$${\alpha }_k$$

$${\beta }_{jk}$$

$${\varphi }_k(x,y)$$

$$=2\frac{\text{precision}\cdot \text{recall}}{\text{precision}+\text{recall}}$$

$$=\frac{\text{\# correctly inferred positions}}{\text{\# inferred positions}}$$

$$=\frac{\text{\# correctly inferred positions}}{\text{\# true positions}}.$$

$$=\exp \left\{-\frac{1}{N}\ln \sum _{\mathbf{Z}}p({\mathbf{X}}_{\text{held-out}},\mathbf{Z}|\mathbf{\Pi },\mathcal{M})\right\}$$

$$=\exp \left\{-\frac{1}{N}\ln p({\mathbf{X}}_{\text{held-out}}|\mathbf{\Pi },\mathcal{M})\right\}$$

Full text

1 Introduction

The alignment of biological sequences (e.g. DNA, RNA and proteins) is one of the most classical and important problems in the field of bioinformatics. Sequence alignment permits the assessment of the functional relationships among biological sequences by quantifying sequence similarity. Because similar nucleotides or amino acids sequences are often functionally related, the development of quantitative evaluations of sequence similarity has been of great interest. This high demand for similarity evaluations has driven the development of a variety of alignment programs [Altschul et al., 1990, Thompson et al., 1994, Frith et al., 2010]. Moreover, sequence alignments are essential for analysing the huge amounts of sequence data produced by high-throughput sequencers in computational tasks such as mapping read sequences onto reference genomes [Li and Homer, 2010, Hamada et al., 2017].

For this alignment task, probabilistic approaches are widely recognised. These probabilistic approaches include Pair Hidden Markov Models (PHMMs) [Durbin et al., 1998], which handle indels and substitutions that occur throughout molecular evolution by using sequentially dependent unobserved hidden states, specifically the match, insertion and deletion states as well as their corresponding probabilistic symbol emissions (cf. Figure 1).

There have been several attempts to construct slightly more complex PHMMs.

Lunter et al. [2008], Bradley and Roberts [2009], Paten et al. [2008] used PHMMs with two insertion and two deletion states, and Cartwright [2009] proposed a general version of PHMMs that employs a zeta power-law model of indel lengths. Additionally, a few generalisations of PHMMs have been proposed, such as Pachter et al. [2002], which introduced generalised PHMMs for DNA–DNA, DNA–cDNA and DNA–protein alignments. However, to the best of our knowledge, no previous study has focused on determining the suitable number of states for representing the biological models that describe sequence evolution or for achieving better alignment accuracy.

Bayesian model selection provides a sophisticated approach for selecting the best model by maximising model evidence. In this, some parameters are marginalised out, and so a preference for simpler models is inherent to the method. When the model prior is uniform, maximising model evidence is equivalent to maximising the posterior probability of model given data, so we can choose the model with the largest posterior by maximising the model evidence.

The well-known difficulty of Bayesian model selection is that the model evidence is analytically intractable in general, including for PHMMs. Markov Chain Monte Carlo (MCMC) [Hastings, 1970] and Variational Inference (VI) [Jordan et al., 1999, Beal, 2003, Blei et al., 2016] enable approximation of the difficult-to-compute model evidence, but both approaches have a drawback: high computational cost. In contrast, the Factorized Asymptotic Bayesian (FAB) algorithm [Fujimaki and Morinaga, 2012, Fujimaki and Hayashi, 2012, Hayashi et al., 2015] is a promising alternative model selection technique based on the Factorized Information Criterion (FIC). One advantage of the FAB algorithm is its simultaneous optimisation of the model structure and the parameters, which makes the FAB algorithm more scalable than VI and MCMC. The advantages are further discussed in section 2.2.

The contributions of this study are summarised as follows.

We developed a novel FIC-based model selection algorithm for PHMMs and demonstrate the reasonably good accuracy in model selection using a synthetic dataset. To the best of our knowledge, this is the first attempt in the literature to apply a model selection method to PHMMs. 2. 2.

The model selection method slightly improved evaluation metrics on the same synthetic dataset. 3. 3.

We conducted experiments on real DNA sequences and found that our method selects a more complex probabilistic structure than the ones that have been traditionally used for pairwise alignment of these species.

2 Methods

PHMMs are a type of probabilistic generative model for sequence alignment [Durbin et al., 1998] with three types of hidden states: a match-type state $M$, an X-insertion-type state $X$ and a Y-insertion-type state $Y$ (Figure 1). The insertion states model the molecular evolution of indels, and the emission probability of the match states characterises the substitution rates.

In this study, we employed the FAB algorithm [Fujimaki and Morinaga, 2012, Hayashi et al., 2015] to select the best model structure for a PHMM. The FAB algorithm is an information criterion-based technique that enables the simultaneous optimisation of both the parameters and the model structure. The properties of the FAB algorithm are explained in 2.2 in more detail. Note that we modified the standard formalisation of PHMM (e.g. Durbin et al. [1998]) because it does not use hidden variables explicitly, which is inappropriate for the FAB algorithm. In this section, we introduce our formalisation of the PHMM with explicit hidden variables (Section 2.1) and then develop a proposed model selection method using the FAB algorithm (Section 2.2).

In the following, we denote the number of match-type states, X-insertion-type states and Y-insertion-type states as $K_{M}$, $K_{X}$ and $K_{Y}$, respectively. Additionally, $K$ represents the total number of hidden states, that is $K=K_{M}+K_{X}+K_{Y}$. Formally, we regard model selection as selecting the number of hidden states $(K_{M},K_{X},K_{Y})$.

2.1 Pairwise Hidden Markov Model (PHMM)

Let observed sequences be $\bm{x}=\{\bm{x}^{n}\}_{n\in[1,N]}$ and $\bm{y}=\{\bm{y}^{n}\}_{n\in[1,N]}$, where $N$ is the number of sequence pairs. The $n$-th sequences are $\bm{x}^{n}=\{x^{n}_{t}\}_{t\in[1,T^{n}_{X}]}$ and $\bm{y}^{n}=\{y^{n}_{u}\}_{u\in[1,T^{n}_{Y}]}$, where $T_{X}^{n}$ and $T_{Y}^{n}$ are the lengths of $\bm{x}^{n}$ and $\bm{y}^{n}$, respectively. We abbreviate all the observed sequences as $\bm{X}=\{\bm{x},\bm{y}\}$. Unlike normal HMMs, PHMMs have hidden variables $\bm{Z}=\{\bm{z}^{n}\}_{n\in[1,N]}$, which are two-dimensional and the $n$-th of which is $\bm{z}^{n}=\{\bm{z}^{n}_{tu}\}_{t\in[0,T^{n}_{X}],u\in[0,T^{n}_{Y}]}$ (Figure 2). Note that these two-dimensional hidden variables are not a common formalisation and include a zero-state, introduced below. The value $\bm{z}^{n}_{t,u}$ corresponds to the hidden state where, for match states, $x^{n}_{t}$ and $y^{n}_{u}$ are matched. For insertion states, $\bm{z}^{n}_{tu}$ represents that $x_{t}$ corresponds to the gap “–” and that the last-used symbol in $\bm{y}^{n}$ is $y^{n}_{u}$ and vice versa for Y-insertion states. The hidden state $\bm{z}_{tu}^{n}=\{z_{t,u,k}^{n}\}_{k\in[1,K]}$ is a 1-of-$K$ representation, but slightly modified to allow a zero-state, where $z^{n}_{t,u,k}$ for all $k$ is zero and does not emit any symbols from that variable (an example is shown in Figure 3). This is because of the unique characteristics of PHMMs (in comparison to conventional HMMs); only a subset of the hidden variables emit symbols, that is only the hidden variables corresponding to aligned positions emit symbols (Figure 3). The set $\bm{\Pi}=\{\bm{\alpha},\bm{\beta},\bm{\phi}\}$ is a parameter set, where $\bm{\alpha},\bm{\beta}$ and $\bm{\phi}$ represent the initial probability, transition probability, and emission probability parameters, respectively. Also, each hidden state $k$ corresponds to one of the state types $\{M,X,Y\}$, which is given by a function $\mathcal{S}$ where $\mathcal{S}(k)\in\{M,X,Y\}$.

Now we can write the complete log-likelihood of PHMM as

[TABLE]

where $\bm{z}^{n}_{\text{in}}$ is a set of hidden variables corresponding to the initial states. The initial hidden variable varies with the type of hidden state because each hidden state’s type uses a different number of original sequences; an $M$-type state uses both $x_{1}$ and $y_{1}$, while $X$- and $Y$-type states each use one of them. For this reason, the initial hidden variable of an $M$-type state is $\bm{z}_{1,1}$, whereas it is $\bm{z}_{1,0}$ and $\bm{z}_{0,1}$ for $X$-type and $Y$-type states, respectively. For the transition probability, we use $pa(\bm{z}^{n}_{tu})=\{z_{t^{\prime}u^{\prime}k}\}_{k\in[1,K]}$ as a set of (previous) hidden variables from which this can transit to $\bm{z}_{t,u}$. We further denote the emission probability as a categorical distribution using new variables $\{\psi^{n}_{tuk}\}$ as

[TABLE]

where $\phi_{k}$ represents the categorical emission probability of the $k$-th hidden state. Note that the $X$- and $Y$-type states emit the gap “–” instead of the normal symbols (e.g. A, T, G or C in the case of DNA alignments). Thus, the dimensionality of the parameter of the emission probability differs with the type of hidden states, namely, $L^{2}-1$ for the match states and $L-1$ for the insertion states, where $L$ is the number of symbols ($L=4$ in the case of DNA sequences). Using this notation, we can rewrite the complete likelihood in an explicit form.

[TABLE]

where $p_{k}(\bm{z}_{tu}|\bm{\beta}_{k})=\prod_{l=1}^{K}p(z_{tul}=1|z_{(t-d_{xl}),(u-d_{yl}),k}=1)^{z_{tul}}$ and $(d_{xk},d_{yk})$ is a transition direction defined as

[TABLE]

Again, it should be noted that the representation in Eq. 3 is essential for a derivation of our model selection algorithm in the following section.

2.2 PHMM model selection algorithm: FAB-PHMM

We formalise the model selection problem for PHMM as a maximisation of the model evidence.

[TABLE]

where the evidence is given by $p(\bm{X}|\mathcal{M})=\int\sum_{\bm{Z}}p(\bm{X},\bm{Z},\bm{\Pi}|\mathcal{M})d\bm{\Pi}$. Note that the model size $\mathcal{M}=\{K_{M},K_{X},K_{Y}\}$ is parameterised by the number of hidden states of each state type. However, the model evidence is difficult to compute; thus, we generally need approximations. In this study, we use FIC as an asymptotically accurate approximation.

FIC has following three appealing properties:

Asymptotic equivalence to marginal likelihood. Although BIC is a widely used and simple information criterion, it lacks theoretical justification because of the non-regularity of the latent variable models [Watanabe, 2009]. PHMM is not an exception to this, so the BIC’s approximation is invalid for PHMMs. Unlike BIC, FIC is consistent with the marginal likelihood for latent variable models. Practically speaking, Fujimaki and Hayashi [2012] empirically showed that BIC-HMM tends to choose overly complicated models, while FIC-HMM chooses optimal models more often. 2. 2.

Simultaneous optimisation of model and parameters. VI is closely related to FIC. Both of them perform similar approximation using variational distribution. One advantage of FIC is it can optimise parameters and models simultaneously. This make FIC-based optimisation computationally more efficient. 3. 3.

Prior free. Unlike VI, FIC does not require prior distributions because it treats priors as $\mathcal{O}(1)$. Thus, FIC is hyper-parameter tuning free and easier to optimise.

In the following, we will start with the derivation of $FIC_{\textit{PHMM}}$ (Section 2.2.1), then take a lower-bound to derive $FICLB$ (Section 2.2.2) for optimisation via expectation maximisation (EM). We iteratively optimise the target function $FICLB$ with respect to a variational distribution $q$ and parameters $\bm{\Pi}$ in the E step (Section 2.2.3) and M step (Section 2.2.4), respectively. The model $\mathcal{M}$ is tuned via model pruning (Section 2.2.5).

2.2.1 FIC

Let $\bm{\Xi}$ be a set of local (component-dependent) parameters and $\bm{\Theta}$ be a set of global (component-independent) parameters. Additionally, we denote all the parameters as $\bm{\Pi}=\{\bm{\Xi},\bm{\Theta}\}$ (in the case of PHMM, local parameters $\bm{\Xi}=\{\bm{\beta},\bm{\phi}\}$ and global parameters $\bm{\Theta}=\{\alpha\}$). Hayashi et al. [2015] have shown that the model evidence can be approximated as an asymptotically accurate information criterion, FIC, which can be expressed as

[TABLE]

where $D_{\bm{\Pi}}$ is the number of free parameters in $\bm{\Pi}$, $\bar{\bm{\Pi}}=\{\bar{\bm{\Xi}},\bar{\bm{\Theta}}\}$ is a maximum joint likelihood estimators (MJLE), $\bm{F}_{\bm{\bar{\Xi}}}$ is the Hessian matrix of $-\ln p(\bm{X},\bm{Z}|\bm{\Pi})$ with respect to $\bm{\bar{\Xi}}$, $q^{**}(\bm{Z})=p(\bm{Z}|\bm{X},\mathcal{M})$ is the marginal posterior and $\mathcal{H}(q^{**})$ is the entropy of $q^{**}$. In FIC, the penalty term is given by the volume of the Fisher information matrix $|F_{\bar{\bm{\Xi}}}|$, which penalises complexity in the model.

Here we derive FIC for PHMM, $FIC_{\textit{PHMM}}$. Since the local parameters $\bm{\Xi}=\{\bm{\beta}_{1},\ldots,\bm{\beta}_{K},\bm{\phi}_{1},\ldots,\bm{\phi}_{K}\}$ do not interact with each other, the Fisher information matrix $F_{\bm{\Xi}}$ is a block diagonal matrix whose blocks are $\{F_{\bm{\beta}_{1}},...,F_{\bm{\beta}_{K}},F_{\bm{\phi}_{1}},\ldots,F_{\bm{\phi}_{K}}\}$, thus $\ln|F_{\Xi}|=\sum_{k}(\ln|F_{\bm{\xi}_{k}}|+\ln|F_{\bm{\beta}_{k}}|)$. Here, using the equation (3), we can write these Fisher information matrices as

[TABLE]

and

[TABLE]

where both $\ln p_{k}(\bm{z}^{n}_{tu}|\bm{\beta}_{k})$ and $\ln p(x^{n}_{t},y^{n}_{u}|\bm{z}^{n}_{tu},\bm{\phi}_{k})$ are $\mathcal{O}(1)$ with respect to the number of samples $N$. Thus, the penalty term is

[TABLE]

The newly introduced symbols $\zeta^{\text{emit}}_{k}$ and $\zeta^{\text{trains}}_{k}$ are effective samples of, respectively, transition and emission probability for the $k$-th latent variable. The values $D_{\bm{\beta}_{k}}$ and $D_{\bm{\phi}_{k}}$ are the dimensionalities of parameters $\bm{\beta}_{k}$ and $\bm{\phi}_{k}$, respectively.

Finaly, ignoring the $\mathcal{O}(1)$ term, we derive FIC for PHMM as

[TABLE]

The penalty terms are now sums of parameter dimensionality weighted by the corresponding effective samples. For example, the dimensionality of the $k$-th emission probability $D_{\phi_{k}}$ is weighted by $\zeta^{\text{emit}}_{k}(\bm{Z})$. When the effective sample of the $k$-th component is small, the penalty term for the $k$-th latent variable also becomes small. In this case, $\bm{Z}$ is degenerate and we can safely prune the $k$-th latent component. This model pruning is further discussed in Section 2.2.5.

2.2.2 FIC Lower-bound

We employ an EM algorithm to optimise the parameters. To make the EM algorithm tractable, we further take the lower bound of $FIC_{\text{PHMM}}$ and derive $FICLB$. We use three approximations to construct the lower bound. (1) Since the MJLEs $\bar{\bm{\Pi}}$ is unavailable in practise, we replace it by the arbitrary parameter $\bm{\Pi}$, which is optimised in the M step. (2) Instead of the marginal posterior $q^{**}$, we use a variational distribution $q$, which is optimised in the E step. (3) We take a lower bound of the negative logarithm as $-\log(\sum_{ntu}z^{n}_{tuk})\geq-L(\sum_{ntu}z^{n}_{tuk},\sum_{ntu}\tilde{q}(z^{n}_{tuk}))$, where $L$ is linear approximation of the logarithm function $L(a,b)=\log b+(a-b)/b$ and $\tilde{q}$ is any distribution over $\bm{Z}$. During the optimisation procedure, $\tilde{q}$ is set to be the variational distribution $q$ of the previous time step. Using these approximations, we now get the lower bound

[TABLE]

Here, we introduced the auxiliary variable $\tilde{\bm{Z}}=\{\tilde{q}(z^{n}_{tuk})\}_{t,u,k,n}$ for simplicity. The full algorithm including model pruning (the model selection mechanism) is explained in Section 2.2.5.

2.2.3 E-step updates

We need to obtain the distribution $q^{*}$ that maximises $FICLB$ (see Algorithm 2.1 for details). This can be done using a modified forward–backward algorithm as follows.

[TABLE]

Using these forward–backward variables, the optimal variational distribution $q^{*}$ is obtained as follows.

[TABLE]

2.2.4 M-step updates

Now, we want to find the $\bm{\Pi}$ that maximises $FICLB$ for fixed $q$ (see Algorithm 2.1 for details). For those parameters, we have the update function

[TABLE]

For calculation of $\beta_{jk}$, out-of-range indexing is treated as zero, that is $q(z^{n}_{d_{xk},d_{yk},k})=0$ if $t-d_{xk}<0$ or $u-d_{yk}<0$.

2.2.5 Pruning degenerated components

In contrast to variational inference, the FAB algorithm enables simultaneous optimisation of a model $\mathcal{M}$ and its parameters $\bm{\Pi}$ via model pruning [Hayashi et al., 2015]. Let us call $\bm{Z}$ degenerated when there exists an equivalent likelihood for a smaller model $\tilde{\mathcal{M}}$, that is $p(\bm{X},\bm{Z}|\bm{\Pi},\mathcal{M})=p(\bm{X},\tilde{\bm{Z}}|\tilde{\bm{\Pi}},\tilde{\mathcal{M}}),$ where $\mathcal{M}>\tilde{\mathcal{M}}$. In such cases, $\mathcal{M}$ is overcomplete so we can transform the model $(\bm{Z},\bm{\Pi})\rightarrow(\tilde{\bm{Z}},\tilde{\bm{\Pi}})$ to obtain a new smaller and equivalent model. This transformation is called model pruning. In the case of FAB-PHMM, we can prune the components $k$ with effective samples $\sum_{ntu}q(z^{n}_{tuk})/N$ that are beneath some threshold $\epsilon$. Starting from a sufficiently large model, we can prune redundant components while optimising parameters. This algorithm is shown in Algorithm 2.1.

We observed that this model pruning algorithm sometimes fails by being captured within poor local optima. In such cases, degenerated components are not pruned. To avoid this problem, we incorporate greedy pruning (Algorithm 2.2). When the algorithm converged, we append the current model $(\bm{\Pi},\mathcal{M},q)$ to model candidates. Then, delete the component with the fewest effective samples (greedy_pruning) and restart the algorithm. After finding the smallest possible model $(K_{M},K_{X},K_{Y})=(1,1,1)$, we choose the model with the largest FIC from among the model candidates.

2.2.6 Computational complexity

For each iteration in Algorithm 2.1, the computational complexity is $\mathcal{O}(N\max(T_{X}^{n})\max(T_{Y}^{n})K^{2})$ for the E and M steps, and $\mathcal{O}(N\max(T_{X}^{n})\max(T_{Y}^{n})K)$ for model pruning. Therefore, the overall complexity for each step is $\mathcal{O}(N\max(T_{X}^{n})\max(T_{Y}^{n})K^{2})$. Note that this complexity is exactly the same as ordinary parameter learning for PHMM using the Baum–Welch algorithm [Durbin et al., 1998].

3 Experiments

We performed three types of experiments to answer the following three questions. First, how accurate is the proposed method in selecting the optimal model (Section 3.1)? Next, how much does the proposed method contribute to alignment accuracy relative to fixed-size models (Section 3.2)? Finally, what kinds of models are selected when we train our proposed method against real DNA data (Section 3.3)?

In this study, because of the high computational cost of parameter learning, we concentrated on short alignments (up to 200 bp). Additionally, we considered only global alignments. Extensions for the alignment of longer sequences and local alignment will be discussed later (in the Discussion section). Moreover, following previous research [Lunter et al., 2008, Bradley and Roberts, 2009], we here concentrate on models with a single match state and multiple insertion states, although our method is potentially applicable to multiple match states (cf. Figure 1 and the Discussion section). For all the experiments, we set the stopping threshold $\eta=10^{-5}$ and pruning threshold $\epsilon=10^{-4}$.

3.1 Model selection capability

We first investigated the model selection capability of the proposed method. We used synthetic data because the true model is not available for real dataset. We defined parameters of PHMMs of different sizes manually and generated alignments from them. The true model size $\mathcal{M}_{\text{true}}=(K^{*}_{M},K^{*}_{X},K^{*}_{Y})$ and their names are shown in Table 1. From the manual models, we generated $N$ alignments of fixed length $100$. After removing the gaps from each alignment, we fed the sequences to Algorithm 2.2 and estimated the optimal model only from the data (pairs of sequences). We then determined if it can recover the true model. We set the initial model size to be $(K_{M},K_{X},K_{Y})=(1,10,10)$. We ran experiments for sample sizes $N=100,200,300,400,500,600,700,800,900,1000$. Parameters for each model can be found in supplementary section S3 and Figures S1–S7.

Table 2 shows the fraction of models that were correctly predicted – prediction is correct when the true model size and predicted model size are exactly the same. As the number of samples $N$ grew, the approximation became more precise. Also, some models tended to require more samples for precise model selection. With sufficient samples, FAB-PHMM successfully recovered models almost perfectly except in some cases (i.e. in experiments (med,$N=1000$) and (small,$N=1000$)). Even in these cases, the inaccurately selected models have smaller FIC values than accurately selected models, and we assumed the algorithm was trapped in a poor local optimum. However, we can easily avoid this problem by using multiple runs and choosing the model with the largest FIC. Indeed, when we picked such a model out of 10 replicates generated with a different random seed, the predictions were always correct when the sample size was greater than or equal to 700 (see bold-face text in Table 2; the bold-face font indicates that the model with the highest FIC value successfully predicted the correct model).

Note that this experimental setting is much simpler than the real setting, where the true model is not in the PHMM class (i.e. the real DNA sequences are not PHMM generated). However, we assume that FAB-PHMM can select the optimal model, in the sense of choosing the closest possible models.

3.2 Alignment accuracy

We also explored the alignment accuracy of the proposed method. Since it is difficult to obtain true genome alignments, we re-used the generated data from the manual models in the model selection experiment. For alignment accuracy assessments, we had five alignment datasets for each different model size. We first trained multiple models on each of those datasets and then performed alignments. For training, we used PHMM (with a fixed model size) of those eight different model sizes with random parameter initializations in addition to FAB-PHMM, which automatically chose the optimal model size. In total, we used nine PHMMs, including one with the true model sizes as well as one FAB-PHMM. In order to avoid poor local optima, we ran five trainings for each setting and selected the model with the best score (FIC for FAB-PHMM and likelihood for PHMM).

As a measure of performance in terms of accuracy, we used the f1 score, which is the harmonic mean of precision and recall:

[TABLE]

For example, when the true alignment is $\genfrac{}{}{0.0pt}{}{x_{1}x_{2}x_{3}x_{4}}{-\ y_{1}y_{2}y_{3}}$ and the inferred alignment is $\genfrac{}{}{0.0pt}{}{x_{1}x_{2}x_{3}x_{4}}{y_{1}\ -y_{2}y_{3}}$, the true and inferred positions are $\{(x_{2},y_{1}),(x_{3},y_{2}),(x_{4},y_{3})\}$ and $\{(x_{1},y_{1}),(x_{3},y_{2}),(x_{4},y_{3})\}$, respectively. The correctly inferred position is simply the intersection of the true positions and inferred positions: $\{(x_{3},y_{2}),(x_{4},y_{3})\}$. In this case, precision is 2/3 and recall is 2/3; thus, the f1 score is 2/3.

We report the result of the number of sequences where $N=1000$ in Table 3. For every dataset, the proposed method performed on par with the true model while all the other fixed-size models performed relatively poorly in some cases, for example for the large model dataset experiment, FAB-PHMM and the large model (i.e. the same model size as the one that produced dataset) performed better than others, whereas the large model performed relatively poorly for smaller datasets.

Although this alignment accuracy measure is widely used (e.g. in Rivas and Eddy [2015]), this approach only considers aligned bases that correspond to match states and ignores all those corresponding to insertion states. For this reason, we also evaluated the insertion counterparts of the f1 score. (See Supplementary Section S1 for further detail.)

In addition to assessing alignment accuracy, we also calculated the perplexity of each trained model in order to show how well the models explain the data. Refer to Supplementary Section S2 for further detail.

3.3 Model selection from real DNA sequences

We used real DNA data to explore the resulting models selected by FAB-PHMM. Because this study concentrates on global alignments of short sequences, proper analyses require homologous DNA sequence pairs. For this reason, we used locally-aligned DNA data produced by Frith and Kawaguchi [2015]111https://zenodo.org/record/17436##.WA3REpOLQYM and multiple-aligned DNA data produced by MULTIZ222http://hgdownload-test.cse.ucsc.edu/goldenPath/hg38/multiz20way/.

3.3.1 LAST dataset

The LAST dataset [Frith and Kawaguchi, 2015] contains pairwise alignments of human sequences to those of four other species (dog, orangutan, mouse and chimpanzee). For our purposes, we selected the alignments with lengths between 100 and 200. Additionally, we removed the alignments with a ”missmap” probability (an alignment ambiguity measure provided with the dataset) of less than $10^{-5}$ because we wanted to use only highly-reliable homologous pairs. This resulted in 1640 to 94,904 remaining alignments for each of the alignments between the human sequences and those of the four species. We extracted overlapping regions from across all four species and cropped alignments to include only those overlapping regions. Then, we randomly sampled 1000 alignments and removed gaps from them in order to input them into FAB-PHMM. We ran 10 trainings with different random seeds (which determine the initial parameters in PHMM) and selected the best model with the highest FIC value. We set the initial model size to be $(K_{M},K_{X},K_{Y})=(1,12,12)$.

Figure 4 shows the selected model size for each species. For species more closely related to humans (i.e. orangutan and chimpanzee), the algorithm selected a similar model and vice versa for species more distantly related to humans (i.e. dog and mouse). Specifically, in the case of chimpanzee and orangutan alignments, the simplest model $(K_{M},K_{X},K_{Y})=(1,1,1)$ was chosen, while many more X- and Y-insertion states were estimated in the case of dog and mouse alignments. Additionally, it was observed that more X-insertion states were predicted than Y-insertion states for dog and mouse alignments.

Figure 5, Figure 6, Supplementary Figure S8 and Supplementary Figure S9 provide the detailed parameters of the selected models from human–chimpanzee, human–mouse, human–orangutan and human–dog sequence alignments, respectively. Overall, all the trained substitution matrixes are almost symmetric. We have the following observations about the human–mouse alignment: two $X$-insertion states, X3 and X4, with similar and relatively high transition probabilities from the match state (M) were predicted, and X3 had higher emission probabilities of A and C while X4 has higher emission probabilities of T and G. In contrast, the human–dog sequence alignment provides the following observations: (1) the self-transition probabilities of X1, X3 and X5 were similar (about 0.8), though the emission probability of X1 had an almost uniform distribution while those of X3 and X5 were skewed with different profiles; (2) the self-transition probability of X2 and X4 were smaller (about 0.4) than those of the others, meaning the states corresponded to shorter gaps; (3) we obtained two long insertion states, Y1 and Y3, and one short insertion state, Y2, while Y2 and Y3 had similar emission profiles.

3.3.2 MULTIZ dataset

We also used the multiz20way dataset, which consists of multiple alignments of 19 different species’ genome assemblies to the human genome. As with the previous experiment, we randomly select 1000 alignments of human sequences to sequences from each of the 19 species where each was restricted to have lengths of 90 to 110 bp. These alignments were then used as a training set in FAB-PHMM. We set the initial model size to be $(K_{M},K_{X},K_{Y})=(1,12,12)$.

The selected model sizes are shown on a phylogenetic tree in Figure 7. In general, we can see species that are more distantly related to humans tend to have alignments with humans that are described by large models, for example the human–bushbaby alignment model has a model size of $(1,9,7)$ and the human–dog alignment model is $(1,9,9)$. However, for species that are more closely related to humans, somewhat random model sizes were selected to describe alignments, for example the human–gorilla alignment model size is $(1,3,3)$ whereas the human–chimpanzee model (a closer species’ alignment model) has a model size of $(1,5,4)$.

4 Discussion

Although our proposed method may be potentially adapted to have multiple match states (cf. Figure 1), we have concentrated on a method with a single match state and multiple ($X$- and $Y$-) insertion states. This is because some previous studies focused on only multiple insertion states, and learning multiple match states in addition to multiple insertion states would lead to more complex models that are not interpretable. (Also, the use of multiple match states incurs more computational cost due to the increase in trained parameters.) Still, it might be interesting to investigate multiple match states for further improvement of alignment accuracy or making novel inferences about sequence evolution, because multiple match states could correspond to the substitution rate of regions.

Our experiments in Section 3.2 show that models trained by our proposed method achieved better alignment accuracy than other models, but the improvement was only marginal. We assume this is because maximising model evidence not always result in maximising alignment accuracy: higher model evidence might contribute to better sequence modelling, but it might not affect alignment accuracy metric. This result is consistent with previous research by Lunter et al. [2008] in which the authors indicated that modifications of insertion states can result in only small improvements in alignment accuracy.

In our analyses using real data (Section 3.3), we utilised limited datasets owing to the high computational cost of our proposed method. Even with these limited datasets, we observed that much more complex models than traditionally used models are selected as optimal. This result implies a possibility that much more complex probabilistic structures exist behind the probabilistic alignments than previously believed.

As well as improving alignment accuracy, the selected model structure may provide interesting insight about biological sequences because the selected probabilistic model structure contains latent information of input data. In other words, the selected model may provide insights into the biological functions of sequences. Similarly, decoded hidden states of PHMMs may reveal (e.g. in DNA alignments) that regions with different hidden states correspond to different functions, such as exon versus intron sequences, non-coding RNAs and regulatory elements. This information can be useful for inferring novel biological insights from sequences.

For the best model selection for real data, however, a comprehensive dataset (including e.g. coding, non-coding and repetitive regions) is required because the trained model structures depend on input data. Indeed, our analyses led to different model selections based on the LAST and MULTIZ datasets. This could be because the homologous pairs taken for the LAST dataset were more similar than those taken for the MULTIZ dataset due to the protocols of generating input sequences (cf. Section 3.3).

In our future work, we will utilise larger and unbiased datasets in order to select more reliable models. The main bottleneck is the high computational cost of our algorithm. To address this, we will attempt to accelerate the training process, such as by parallelization of the algorithms, stochastic optimisation [Robbins and Monro, 1951, Hoffman et al., 2013, Liu et al., 2015] and seed-extension heuristics in forward-and-backward algorithms [Hamada et al., 2017]

In this study, we focused on genomic DNA sequences, but our method is applicable to RNA or protein sequences as well. The number of characters in protein sequences greatly exceeds those in DNA sequences, which would lead to more complex models that may provide interesting biological insights. These applications will be included in our future research as well.

5 Conclusion

In this study, we proposed a novel method to develop PHMMs based on FIC and demonstrated the model selection capability of the proposed model using a synthetic dataset. We believe this is the first study that focuses on model selection of PHMM. On the same synthetic dataset, we observed slight improvement of evaluation metrics of sequence alignments. Additionally, we conducted experiments on real DNA sequences and found that they are best handled with a more complex probabilistic structure than the ones that have been traditionally used for pairwise alignment of these species. This result implies a possibility that more complex probabilistic structures exist behind probabilistic alignments than previously believed.

Acknowledgement

The authors are grateful to Kohei Hayashi (AIST), Tsukasa Fukunaga (Waseda University) and members of the Hamada Laboratory at Waseda University. Computation for this study was partially performed on the NIG supercomputer at ROIS National Institute of Genetics.

Funding

This work was supported by MEXT KAKENHI Grant Numbers JP24680031, JP16H05879 and JP25240044 to MH, in part.

Appendix S1 Alignment accuracy evaluation for insertion states

In addition to measuring accuracy according to aligned positions (explained in Section 3.2), we also evaluated insertion accuracy. We define this measure as precision/recall/f1-score of inserted. Using the same example as used in Section 3.2, when the true alignment is $\genfrac{}{}{0.0pt}{}{x_{1}x_{2}x_{3}x_{4}}{-\ y_{1}y_{2}y_{3}}$ and the inferred alignment is $\genfrac{}{}{0.0pt}{}{x_{1}x_{2}x_{3}x_{4}}{y_{1}\ -y_{2}y_{3}}$, true insertions and inferred insertions are $\{(x_{1},-)\}$ and $\{(x_{2},-)\}$. In this case, all of precision/recall/f1 are 0. As shown in Table S1, this measure is stricter than traditional measures (Table 3) because predicting matched positions is easier than predicting insertions.

Appendix S2 Perplexity of trained models

Perplexity is a measure of how accurately a probabilistic model predicts a sample (in this case a sample properly corresponding to an alignment) and is defined as

[TABLE]

where $\bm{X}_{\text{held-out}}$ represents the reserved sequences (i.e. sequences not used for training) and $N$ is the number of sequences (500 sequences in this experiment).

Because perplexity can be very large, we report log perplexity instead (Supplementary Table S2). The proposed method outperformed smaller trained models when simulated models were large. For example, when the simulated model was huge, our proposed method achieved a log perplexity of 177.0064 but the med model’s log perplexity is only 179.7014. Because small or med models are usually considered for alignment tasks, our proposed model has advantages for alignment modelling over existing methods. It is unexpected but interesting that larger models (such as huge and large) achieve accuracy that is comparable to that of the proposed method, which indicates that the larger models do not result in overfitting even when the data are generated by simpler models.

Appendix S3 Parameters for artificial models

Here, we show initial/transition/emission probabilities of each artificial model in Figure S1 (small), S2 (med), S3 (large), S4 (huge), S5 (imb), S6 (imb_large) and S7 (imb_huge).

Appendix S4 Parameter visualisation of the model inferred from the LAST dataset