- Proceedings
- Open Access

# Statistical tests to identify appropriate types of nucleotide sequence recoding in molecular phylogenetics

- Victor A Vera-Ruiz
^{1}, - Kwok W Lau
^{2, 3}, - John Robinson
^{1}and - Lars S Jermiin
^{2, 4}Email author

**15 (Suppl 2)**:S8

https://doi.org/10.1186/1471-2105-15-S2-S8

© Vera-Ruiz et al.; licensee BioMed Central Ltd. 2014

**Published:**31 January 2014

## Abstract

### Background

Under a Markov model of evolution, recoding, or lumping, of the four nucleotides into fewer groups may permit analysis under simpler conditions but may unfortunately yield misleading results unless the evolutionary process of the recoded groups remains Markovian. If a Markov process is lumpable, then the evolutionary process of the recoded groups is Markovian.

### Results

We consider stationary, reversible, and homogeneous Markov processes on two taxa and compare three tests for lumpability: one using an *ad hoc* test statistic, which is based on an index that is evaluated using a bootstrap approximation of its distribution; one that is based on a test proposed specifically for Markov chains; and one using a likelihood-ratio test. We show that the likelihood-ratio test is more powerful than the index test, which is more powerful than that based on the Markov chain test statistic. We also show that for stationary processes on binary trees with more than two taxa, the tests can be applied to all pairs. Finally, we show that if the process is lumpable, then estimates obtained under the recoded model agree with estimates obtained under the original model, whereas, if the process is not lumpable, then these estimates can differ substantially. We apply the new likelihood-ratio test for lumpability to two primate data sets, one with a mitochondrial origin and one with a nuclear origin.

### Conclusions

Recoding may result in biased phylogenetic estimates because the original evolutionary process is not lumpable. Accordingly, testing for lumpability should be done prior to phylogenetic analysis of recoded data.

## Keywords

- Phylogeny
- Markov model
- stationarity
- homogeneity,reversibility
- recoding
- lumping
- nucleotides
- primates

## Introduction

When nucleotides intentionally are recoded to a 3- or 2-state alphabet in order to focus on a subset of the possible types of substitutions (e.g., transversions [1–3]) or reduce compositional heterogeneity [4], it is no longer appropriate to use model-based phylogenetic methods that rely solely on time-reversible, 4-state Markov models. Instead, one needs to use a 3- or 2-state Markov model to approximate the evolutionary processes for the recoded sequence data. This requirement was first realised by Phillips and Penny [5], who used a time-reversible 2-state Markov model [6] to analyse *RY*-recoded nucleotide sequences, and Gibson et al. [7], who developed a time-reversible 3-state Markov model to analyse *Y*-recoded nucleotide sequences. Before these studies, other investigators had used *RY*-recoded nucleotide sequences to infer the evolutionary relationships among mammals [1–3] and among bacteria [4].

Recoding of nucleotides and/or amino acids has been used repeatedly in recent phylogenetic studies [8–31]. However, the mathematical principles underpinning the recoding of nucleotides or amino acids have not yet been adequately examined. For example, it is not yet known whether the Markovian property is maintained after recoding and how this should be tested [32]. Without this knowledge, we may run the risk of using a promising procedure in a manner that turns out to be inappropriate for the data.

In this paper, we take a first step by considering tests for lumpability in a Markov model of evolution for pairs of homologous nucleotide sequences (we are aware of only one paper in the phylogenetic literature where the term lumpability is used [33], but there it was used in a different context). We only consider nucleotides but believe our tests could be generalized to encompass amino acids as well. We then illustrate the performance of our tests for lumpability using simulated and real data, and show that recoding of nucleotides should be used with caution when analysing DNA phylogenetically.

## Methods

### The theoretical basis for recoding nucleotides

*q <*4. Then $\phantom{\rule{0.5em}{0ex}}\mathcal{S}$ is reduced by grouping, or lumping, some of the original states (i.e.,

*A*,

*C*,

*G*, and

*T*) into one or two new states (i.e.,

*R*,

*Y*,

*S*,

*W*,

*M*,

*K*,

*B*,

*D*,

*H*, and

*V*)--in molecular phylogenetics, this procedure has been called

*recoding*[34]. Table 1 presents the 13 possible recoding schemes and partitions of ${\mathcal{S}}^{\prime}$ using notation established by the NC-IUB [35]. The 13 recoding schemes fall into three major grouping categories, as shown in Table 1.

The 13 ways of reducing a 4-letter state space ($\phantom{\rule{0.5em}{0ex}}\mathcal{S}$) to a 3- or 2-letter state space ${\mathcal{S}}^{\prime}$.

Nucleotide-grouping Subsets | Recoding notation | Resulting ${\mathcal{S}}^{\prime}$ | Major grouping category |
---|---|---|---|

{{ |
| { | {2 : 1 : 1} |

{ |
| { | |

{ |
| { | |

{ |
| { | |

{ |
| { | |

{ |
| { | |

{ |
| { | {3 : 1} |

{ |
| { | |

{ |
| { | |

{ |
| { | |

{{ |
| { | {2 : 2} |

{{ |
| { | |

{{ |
| { |

### The evolutionary process for two homologous nucleotide sequences

*n*independently evolving sites, which have diverged under Markovian conditions from their common ancestor on a rooted, 2-tipped tree. Let

*π*

_{0}denote the initial probability vector of the nucleotide frequencies, such that ${\pi}_{0}^{T}=\left({\pi}_{{0}_{1}},{\pi}_{{0}_{2}},{\pi}_{{0}_{3}},{\pi}_{{0}_{4}}\right)$, where, for convenience of notation, we will use the subscripts 1, 2, 3, 4 to denote

*A*,

*G*,

*C*,

*T*. Over each edge of this tree, there is a substitution process,

*X*(

*t*) and

*Y*(

*t*), respectively, described by the transition probabilities

*f*

_{ ij }(

*t*) denote the theoretical joint probability of a site being in state

*i*in A and state

*j*in B at time

*t*:

**F**(

*t*) = {

*f*

_{ ij }(

*t*)} denote the joint probability matrix and let

**P**

^{ X }(

*t*) and

**P**

^{ Y }(

*t*) denote the transition probability matrices of

*X*(

*t*) and

*Y*(

*t*). In practice, the two matrices cannot be identified from

**F**(

*t*) without some assumptions about the evolutionary processes of the sequences A and B. We assume that the processes are globally stationary, reversible, and homogeneous (SRH) (for definitions, see [36–39], and, in more detail, [40]). Given these assumptions, take

**P**

^{ X }(

*t*) =

**P**

^{ Y }(

*t*) =

**P**(

*t*) and, if

*π*

_{ X }and

*π*

_{ Y }denote the equilibrium probability distributions of the processes

*X*(

*t*) and

*Y*(

*t*), take

*π*

_{ X }=

*π*

_{ Y }=

*π*

_{0}=

*π*and write Π =

*diag*(

*π*). Then, from (1), we get

**F**(

*t*) =

**P**(

*t*)

^{ T }Π

**P**(

*t*). The transition probability matrix can be expressed by an instantaneous rate matrix

**R**, such that

**P**(

*t*) =

*e*

^{ R t }, where

*R*

_{ ij }

*≥*0 for

*i ≠ j*, ${R}_{ii}=-{\sum}_{i\ne j}^{4}{R}_{ij},$ and

*π*

^{ T }

**R**=

**0**

^{ T }, where

*π*

^{ T }is the equilibrium distribution of

**R**[40]. Furthermore, the instantaneous rate matrix can take the form

**R**=

**S**Π, where

**S**is a symmetric matrix with

*s*

_{ ij }

*≥*0 for

*i ≠*

*j*, and ${s}_{ii}=-{\sum}_{i\ne j}^{4}{s}_{ij}{\pi}_{j}\mathsf{\text{/}}{\pi}_{i}$[40]. The matrices

**R**and

**P**(

*t*) can be written in terms of the eigenvector decomposition of Π

^{1/ 2}

**S**Π

^{1/ 2}=

**L**Λ

**L**

^{ T }. In other words

**Λ**is a diagonal matrix with columns containing the eigenvalues of

**Π**

^{1/ 2}

**SΠ**

^{1/ 2}and

**L**is a matrix with columns containing its right eigenvectors. The joint probability matrix is then symmetric and

Note that under these assumptions, there are only nine free parameters to be estimated: six free parameters for the off-diagonal elements of **S** (define *s*^{
T
} = (*s*_{12}, *s*_{13}, *s*_{14}, *s*_{23}, *s*_{24}, *s*_{34})) and three free parameters for *π* (because ${\pi}_{4}=1-{\sum}_{i=1}^{3}{\pi}_{i}$). The time *t* can be fixed at 1 since modifying it is equivalent to modifying the *s*-parameters.

**N**denote the 4

*×*4 divergence matrix for A and B, such that

**N**= {

*n*

_{ ij }}, where

*n*

_{ ij }represents the number of homologous sites that are in state

*i*in A and state

*j*in B. Under the model, the vector of elements of

**N**has a multinomial distribution with parameters

*n*and

**F**(

*t*); its expected value is thus

*E*(

**N**) =

*n*

**F**(

*t*). Because the parameters

*s*and

*π*are in a one-to-one relation with the elements of

**F**, the maximum-likelihood estimates of

*s*and

*π*can be obtained from the eigenvector decomposition of $\widehat{F}\left(t\right)=\frac{1}{2n}\left(N+{N}^{T}\right)$, then

where $\widehat{L}$ and $\widehat{\text{\Lambda}}$ are obtained from ${\widehat{\mathbf{\text{\Pi}}}}^{-1/2}\widehat{F}\left(1\right){\widehat{\mathbf{\text{\Pi}}}}^{-1/2}=\widehat{L}{e}^{2\widehat{\mathbf{\text{\Lambda}}}}{\widehat{L}}^{T}$.

### Lumpable Markov chains

*q <*4, such that a

*lumped process*, ${X}^{\prime}\left(t\right)$, with a smaller number of states, is generated with transition probabilities

and initial probabilities ${\pi}^{\prime}=P\left[{X}^{\prime}\left(0\right)=k\right]=P\left[X\left(0\right)\in {\mathcal{S}}_{k}\right]$.

*π*, the lumped process, defined in (7), is a Markov chain whose transition probabilities do not depend on the choice of

*π*. A necessary and sufficient condition for

*X'*(

*t*) to be lumpable with respect to a partition ${\mathcal{S}}^{\prime}$ is that for every pair of subsets, ${\mathcal{S}}_{k}$ and ${\mathcal{S}}_{l}$, ${\sum}_{j\in {S}_{l}}{P}_{ij}\left(t\right)$, has the same value for every state

*i*in ${\mathcal{S}}_{k}$[41, 42]. Accordingly, if

*X'*(

*t*) is lumpable, then the transition probabilities for

*X'*(

*t*) for any given pair of subsets in ${\mathcal{S}}^{\prime}$ are

**P**'(

*t*) can be expressed as a matrix function of

**P**(

*t*) as follows:

**V**is a 4

*× q*matrix, where

*q*is the number of states in the lumped process, such that the

*l*-th column of

**V**is a vector with 1's in the components corresponding to states in

*S*

_{ l }and 0's otherwise, and

*q ×*4 matrix whose

*k*-th row is a probability vector with non-zero elements corresponding to the states in ${\mathcal{S}}_{k}$. A useful necessary and sufficient condition for lumpability [41, 42] is

Conditions required for a 4-state Markovian process to be lumpable (in terms of *s* and *π*), and transformations to obtain $\tilde{\pi}$ and $\tilde{s}$ such that the lumpability holds.

${\mathcal{S}}^{\prime}$ | Lumpability conditions | ($\tilde{\pi}$,$\tilde{s}$) |
---|---|---|

{{ |
| ${\tilde{s}}_{\mathsf{\text{13}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{23}}}=\left({\u015d}_{\mathsf{\text{13}}}+{\u015d}_{\mathsf{\text{23}}}\right)/\mathsf{\text{2}}$ ${\tilde{s}}_{\mathsf{\text{14}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{34}}}=\left({\u015d}_{\mathsf{\text{14}}}+{\u015d}_{\mathsf{\text{34}}}\right)/\mathsf{\text{2}}$ |

{ |
| ${\tilde{s}}_{\mathsf{\text{12}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{14}=\left({\u015d}_{\mathsf{\text{12}}}+{\u015d}_{\mathsf{\text{14}}}\right)/\mathsf{\text{2}}$ ${\tilde{s}}_{\mathsf{\text{23}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{34}}}=\left({\u015d}_{\mathsf{\text{23}}}+{\u015d}_{\mathsf{\text{34}}}\right)/\mathsf{\text{2}}$ |

{ |
| ${\tilde{s}}_{\mathsf{\text{12}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{13}}}=\left({\u015d}_{\mathsf{\text{12}}}+{\u015d}_{\mathsf{\text{13}}}\right)/\mathsf{\text{2}}$ ${\tilde{s}}_{\mathsf{\text{24}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{34}}}=\left({\u015d}_{\mathsf{\text{24}}}+{\u015d}_{\mathsf{\text{34}}}\right)/\mathsf{\text{2}}$ |

{ |
| ${\tilde{s}}_{\mathsf{\text{12}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{24}}}=\left({\u015d}_{\mathsf{\text{12}}}+{\u015d}_{\mathsf{\text{24}}}\right)/\mathsf{\text{2}}$ ${\tilde{s}}_{\mathsf{\text{13}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{34}}}=\left({\u015d}_{\mathsf{\text{13}}}+{\u015d}_{\mathsf{\text{34}}}\right)/\mathsf{\text{2}}$ |

{{ |
| ${\tilde{s}}_{\mathsf{\text{13}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{23}}}=\left({\u015d}_{\mathsf{\text{13}}}+{\u015d}_{\mathsf{\text{23}}}\right)/\mathsf{\text{2}}$ ${\tilde{s}}_{\mathsf{\text{14}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{24}}}=\left({\u015d}_{\mathsf{\text{14}}}+{\u015d}_{\mathsf{\text{24}}}\right)/\mathsf{\text{2}}$ |

{ |
| ${\tilde{s}}_{\mathsf{\text{13}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{14}}}=\left({\u015d}_{\mathsf{\text{13}}}+{\u015d}_{\mathsf{\text{14}}}\right)/\mathsf{\text{2}}$ ${\tilde{s}}_{\mathsf{\text{23}}}={\tilde{s}}_{\mathsf{\text{24}}}=\left({\u015d}_{\mathsf{\text{23}}}+{\u015d}_{\mathsf{\text{24}}}\right)/\mathsf{\text{2}}$ |

{ |
| ${\tilde{s}}_{\mathsf{\text{12}}}={\tilde{s}}_{\mathsf{\text{13}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{14}}}=\left({\u015d}_{\mathsf{\text{12}}}+{\u015d}_{\mathsf{\text{13}}}+{\u015d}_{\mathsf{\text{14}}}\right)/\mathsf{\text{3}}$ |

{ |
| ${\tilde{s}}_{\mathsf{\text{12}}}={\tilde{s}}_{\mathsf{\text{23}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{24}}}=\left({\u015d}_{\mathsf{\text{12}}}+{\u015d}_{\mathsf{\text{23}}}+{\u015d}_{\mathsf{\text{24}}}\right)/\mathsf{\text{3}}$ |

{ |
| ${\tilde{s}}_{\mathsf{\text{13}}}={\tilde{s}}_{\mathsf{\text{23}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{34}}}=\left({\u015d}_{\mathsf{\text{13}}}+{\u015d}_{\mathsf{\text{23}}}+{\u015d}_{\mathsf{\text{34}}}\right)/\mathsf{\text{3}}$ |

{ |
| ${\tilde{s}}_{\mathsf{\text{14}}}={\tilde{s}}_{\mathsf{\text{24}}}=\phantom{\rule{0.5em}{0ex}}{\tilde{s}}_{\mathsf{\text{34}}}=\left({\u015d}_{\mathsf{\text{14}}}+{\u015d}_{\mathsf{\text{24}}}+{\u015d}_{\mathsf{\text{34}}}\right)/\mathsf{\text{3}}$ |

{{ |
| ${\tilde{s}}_{23}=\frac{{\tilde{s}}_{12}\left({\widehat{\pi}}_{2}{\widehat{\pi}}_{3}-{\widehat{\pi}}_{1}{\widehat{\pi}}_{4}\right)\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}{\tilde{s}}_{14}{\widehat{\pi}}_{4}\left({\widehat{\pi}}_{1}+{\widehat{\pi}}_{3}\right)\phantom{\rule{0.3em}{0ex}}}{{\widehat{\pi}}_{3}\left({\widehat{\pi}}_{2}+{\widehat{\pi}}_{4}\right)\phantom{\rule{0.3em}{0ex}}}$ ${\tilde{s}}_{34}=\frac{{\tilde{s}}_{12}{\widehat{\pi}}_{1}+{\tilde{s}}_{23}{\widehat{\pi}}_{3}-{\tilde{s}}_{14}{\widehat{\pi}}_{1}}{{\widehat{\pi}}_{3}}$ |

{{ |
| ${\tilde{s}}_{13}=\frac{{\tilde{s}}_{12}\left({\widehat{\pi}}_{1}{\widehat{\pi}}_{3}-{\widehat{\pi}}_{2}{\widehat{\pi}}_{4}\right)\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}{\tilde{s}}_{24}{\widehat{\pi}}_{4}\left({\widehat{\pi}}_{2}+{\widehat{\pi}}_{3}\right)\phantom{\rule{0.3em}{0ex}}}{{\widehat{\pi}}_{3}\left({\widehat{\pi}}_{1}+{\widehat{\pi}}_{4}\right)\phantom{\rule{0.3em}{0ex}}}$ ${\tilde{s}}_{34}=\frac{{\tilde{s}}_{12}{\widehat{\pi}}_{2}+{\tilde{s}}_{13}{\widehat{\pi}}_{3}-{\tilde{s}}_{24}{\widehat{\pi}}_{2}}{{\widehat{\pi}}_{3}}$ |

{{ |
| ${\tilde{s}}_{23}=\frac{{\tilde{s}}_{13}\left({\widehat{\pi}}_{2}{\widehat{\pi}}_{3}-{\widehat{\pi}}_{1}{\widehat{\pi}}_{4}\right)\phantom{\rule{0.3em}{0ex}}+\phantom{\rule{0.3em}{0ex}}{\tilde{s}}_{14}{\widehat{\pi}}_{4}\left({\widehat{\pi}}_{1}+{\widehat{\pi}}_{2}\right)\phantom{\rule{0.3em}{0ex}}}{{\widehat{\pi}}_{2}\left({\widehat{\pi}}_{3}+{\widehat{\pi}}_{4}\right)\phantom{\rule{0.3em}{0ex}}}$ ${\tilde{s}}_{24}=\frac{{\tilde{s}}_{13}{\widehat{\pi}}_{1}+{\tilde{s}}_{23}{\widehat{\pi}}_{2}-{\tilde{s}}_{14}{\widehat{\pi}}_{1}}{{\widehat{\pi}}_{2}}$ |

We note that under certain conditions, such as those considered by the JC [43] and F81 [44] models, all recodings are lumpable. Conditions under which recoding of nucleotides are possible for the K2P [45] and HKY [46] models are given in Table 3.2 of [32].

### Tests for lumpability

We consider three possible tests: An *ad hoc* test based on a parametric bootstrap for an index of departure from the lumpability condition [32]; a test based on a test for lumpability in Markov chains [47]; and a likelihood-ratio test.

#### Index test

*η ≥*0, with

*η*being 0 only under lumpable Markovian processes. Then, the hypothesis that the Markov process is lumpable is equivalent to the hypothesis

*H*

_{0}:

*η*= 0. From the observed divergence matrix,

**N**of two homologous sequences, assuming a SRH Markovian model of evolution, an estimate $\widehat{\eta}$ can be used as a test statistic for

*H*

_{0}, where

and $\widehat{\mathbf{M}}=VU\widehat{P}\left(1\right)V-\widehat{P}\left(1\right)V$ for $\widehat{P}\left(1\right)=\mathsf{\text{exp}}\left(\widehat{S}\widehat{\text{\Pi}}\right)$.

*B*matrices can be generated by simulation under conditions of lumpability, where we take ${N}_{b}^{*}$, with

*b*∈ {1,

*⋯*,

*B*}, to be independent and multinomial with parameters

*n*and $\tilde{\mathbf{F}}\left(1\right)$. From each of these simulated samples, we calculate ${F}_{b}^{*}\mathsf{\text{(1)}}=\left({N}_{b}^{*}+{N}_{b}^{*T}\right)\mathsf{\text{/2}}$, ${\text{\Pi}}_{b}^{*}$ and ${S}_{b}^{*}$ from ${F}_{b}^{*}$, as in (5) and (6), and then ${P}_{b}^{*}\left(1\right)=\mathsf{\text{exp}}\left({S}_{b}^{*}{\text{\Pi}}_{b}^{*}\right),{M}_{b}^{*}={VUP}_{b}^{*}\left(1\right)V-{P}_{b}^{*}\left(1\right)V$, and

The true *P*-value is then the probability that we obtain a value as large as or larger than the observed $\widehat{\eta}$, so a bootstrap approximation to this *P*-value is the proportion of ${\eta}_{i}^{*},...,{\eta}_{B}^{*}$ exceeding $\widehat{\eta}$.

#### Markov chain test

*χ*

^{2}test to determine whether a Markov chain is lumpable with respect to a partition ${\mathcal{S}}^{\prime}$ is available [47]. The test is based on the comparison of observed transition frequencies to their respective theoretical counterparts under the null hypothesis that the chain is lumpable. The approach does not make any assumption about reversibility or stationarity of the process. The authors used a matrix of transition counts, {

*n*

_{ ij }}, to estimate the transition probabilities

*p*

_{ ij }, where

*n*

_{ ij }represents the number of transitions into state

*j*from state

*i*in one step, so the number of steps in the Markov chain is

*n*

_{ •• }, where the subscript

*•*indicates summation. Now, if we start from our divergence matrix

**N**, where

*n*

_{ ij }represents the number of sites that are in state

*i*for sequence A and state

*j*in sequence B, and the SRH assumptions are kept, either A or B can be assumed to be the original sequence at time 0, whereas the other one can be assumed to be the observed sequence at time 2 (since we took the edge lengths to be 1). Take A as the ancestral sequence, then the divergence matrix

**N**has the same properties as a transition count matrix, and we can proceed as described in [47]. A transition probability from

*i*to ${\mathcal{S}}_{l}$ is

*l*= 1,

*..*.,

*q*. From the definition of a lumpable process (7), if the Markovian process is lumpable with respect to ${\mathcal{S}}^{\prime}$, then

*γ*

_{ k }is the number of states that are part of the subset ${\mathcal{S}}_{k}$, and

*k*= 1,

*..*.,

*q*. Therefore, ${g}_{il}={P}_{kl}^{\prime}\left(\mathsf{\text{2}}\right)$ if the process is lumpable and the null hypothesis of lumpability can be expressed as

*H*

_{0}: ${g}_{il}={P}_{kl}^{\prime}\left(\mathsf{\text{2}}\right)$ for all

*i*∈ ${\mathcal{S}}_{k}$. Given the divergence matrix,

**N**, estimates

*g*

_{ il }and ${P}_{kl}^{\prime}\left(\mathsf{\text{2}}\right)$ are

and showed (by pointing out that *o*_{
il
} *− e*_{
il
} are a stack of *q* tables of size 4 *× γ*_{
l
} of mean-corrected multinomials with row and column sums equal to zero) that the test statistic is distributed under *H*_{0} as a *χ*^{2} variable with $\left(q-\mathsf{\text{1}}\right){\sum}_{k=1}^{q}\left({\gamma}_{k}-\mathsf{\text{1}}\right)$ degrees of freedom, if all cells are non-zero. In the case considered here, the degrees of freedom for any of the recoding schemes is 2.

#### Likelihood-ratio test

*n*

_{ ij }} =

**N**, the observed divergence matrix,

**F**

_{(π,s)}(1) = exp(

**SΠ**) and {

*f*

_{ ij }(

*π*,

*s*)} =

**F**

_{(π,s)}(1). These matrices are obtained as shown in (5) and (6). We also want to estimate (

*π*,

*s*) under the constraints imposed by the null hypothesis of lumpablity,

*H*

_{0}. The constraints are given in the second column of Table 2. Then we can define the constrained estimates ($\tilde{\pi}$, $\tilde{s}$) to satisfy

**A**, such that

*y*is the response constraint vector, defined such that two values of

*y*are zero corresponding to the two constraints. The matrix

**A**will, in the case of partitions into two groups of two, contain

*π*, so to emphasise this possible dependence, write

*y*=

*g*(

*s|π*). Also write

*s*=

**A**

^{−1}

*y*=

*g*

^{−1}(

*y|π*). Then

The optimization process is done in two steps: the values of *s*, if dependent of *π*, are optimized given the original *π* set, then the *π* vector is optimized given the optimized values of *s*. This process is repeated until convergence is achieved.

*LR*, can be calculated with

Under the null hypothesis of lumpability, 2 *× LR* is distributed as a *χ*^{2} variable with 2 degrees of freedom.

## Results

### Assessment of accuracy

The joint probability distribution was calculated by the steps given in (2), (3), and (4); then, assuming a nucleotide sequence of length *n* = 1500, 5000 divergence matrices were calculated by Monte Carlo simulations assuming that **N**_{
i
} is multinomial with parameters (*n*, **F**(1)) for *i* = 1, *. .* . , 5000.

*PP*plot displaying the distribution of observed

*P*-values, obtained from each test, plotted against the expected

*P*-values, obtained from the uniform distribution. The linear relationship between these two sets of

*P*-values (Figure 1) confirms the accuracy of the tests.

### Comparisons of power

The power of each test was compared for each recoding scheme under non-lumpable conditions. To do this, we used *π*^{
T
} = (0.1, 0.2, 0.3, 0.4) and values of *s* that yield increasing values of *η*, as in (9), generated 3000 divergence matrices using Monte Carlo simulation, and then calculated the three test statistics and their corresponding *P*-values, using the procedures explained above, for each value of *η*. The power at the 5% level, is then equal to the proportion of observed *P*-values less than 0.05.

*RY*recoding--similar power curves were obtained for the other 12 recoding schemes (results not shown). All of these results indicate that the likelihood-ratio test is the most powerful of the tests considered, followed by the Index test and, finally, by the Markov Chain test.

### Cases with more than 2 homologous sequences

For general cases involving more than 2 homologous sequences, we can test for lumpability in all pairs of sequences by using the methods described above under the assumption that the evolutionary process is SRH over the whole tree (i.e., the process is globally SRH). For example, in the case of an alignment with seven sequences, there will be 21 *P*-values. A *PP*-plot with these *P*-values should yield a straight line when the data are lumpable, and deviations from this expectation when the processes are not lumpable. However, the observed *P*-values are not independent, so we need to show that this condition is not cause for concern for the dots in a *PP*-plot to be on a straight line. We give a simplified argument taken from [48]. Consider a set of observed *P*-values *P*_{1}, . . ., *P*_{
n
}, which, if all the null hypotheses are true, are identically distributed as uniform random variables on (0, 1). Let *p* be any value between 0 and 1, let *I*(*P*_{
j
} *< p*) take the value 1 if *P*_{
j
} *< p* and 0 otherwise, and let ${N}_{p}={\sum}_{j=1}^{n}I\left({P}_{j}<p\right)$ be the number of observed *P*-values less than *p*. Then *E*(*I*(*P*_{
j
} *< p*)) = *P* (*P*_{
j
} *< p*) = *p* and so the expected number of *P*-values less than *p* is *E*(*N*_{
p
}) = *np*. This implies that the plot of the observed *P*-values will lie approximately on a straight line. The dependence will cause some clustering of the observed *P*-values but the *PP*-plot will remain useful in indicating whether there is evidence against some of the hypotheses.

*PP*-plots shown in Figure 3 were obtained from alignments of nucleotides generated under lumpable or non-lumpable conditions, with respect to

*RY*recoding, on the tree shown in Figure 4 before being analysed using the likelihood-ratio test. From these two plots, it is clear that the test is able to identify cases where sequences have evolved under non-lumpable conditions.

### The effect of non-lumpability on phylogenetic estimates

*RY*), then we can obtain

**F**

_{ q }(

*t*) =

**V**

^{ T }P(

*t*)

**V**and from this, using (2, 3) and (4), we can further obtain the transition matrix of the process

*X'*(

*t*),

where **Π**_{
q
} = **V**^{
T
}**ΠV** . If the process is not lumpable with respect to that recoding scheme, then *X'*(*t*) is not a Markov process and, although we can calculate **P**_{
q
}(*t*), it is *not* the transition matrix of the process *X'*(*t*).

In either case, the matrix **P** *'*(*t*) = **UP**(*t*)**V** can be defined and, if the process *X*(*t*) is lumpable, then *P'*(*t*) = *P*_{
q
}(*t*). On the other hand, if *X*(*t*) is not lumpable, the elements of **P** *'*(*t*) are still given as ${\mathbf{P}}_{kl}^{\prime}\mathsf{\text{(}}t\mathsf{\text{)}}=\mathbf{P}\left({X}^{\prime}\mathsf{\text{(}}t\mathsf{\text{)}}=l\mathsf{\text{|}}{X}^{\prime}\mathsf{\text{(}}0\mathsf{\text{)}}\phantom{\rule{2.77695pt}{0ex}}\mathsf{\text{=}}\phantom{\rule{2.77695pt}{0ex}}k\right)$, but **P** *'*(*t*) is no longer a transition matrix of *X'*(*t*). Conveniently, we can compare **P** *'*(1), the true conditional probability matrix at *t* = 1, with **P**_{
q
} (1), the false transition matrix at *t* = 1, thus allowing us to examine the effect of non-lumpability.

*π*

^{ T }= (0.2, 0.3, 0.2, 0.3) and

*s*

^{ T }= (0.2, 0.1, 0.3, 0.3, 1.0, 0.2). As can be seen from the

*π*- and

*s*-vectors, the lumpable condition is met for

*RY*-recoding but not for

*KM*-recoding. Figures 5b and 5c display the corresponding tree with the edge lengths adjusted according to, respectively, the

*RY*- and

*KM*-recoding schemes.

*RY*-recoded data is shorter than the corresponding edge in the tree obtained from the original data. However, because the original process was lumpable with respect to

*RY*-recoding, the relative length of each edge in the two trees is the same, the difference being equal to a scale factor

Every edge in the tree obtained from the *KM*-recoded data is also shorter than the corresponding edge in the tree obtained from the original data, but the relative length of each edge in the two trees differ, the reason being that the process generating the original data was not lumpable with respect to *KM*-recoding.

Figures 5d, 5e, and 5f show the corresponding results for the data generated by Monte Carlo simulation. The three trees display the same characteristics as those shown in Figures 5a, 5b, and 5c, while also showing some variation in the edge lengths that is due to the random nature of the data and the finite sample size. Hence, although recoding of nucleotides might be useful for a variety of reasons, using recoded data, without having tested for lumpability first, might lead to biased phylogenetic estimates.

### Example 1 -- Primate mitochondrial DNA

*PP*-plots from tests for lumpability for

*RY*recoding, indicating non-lumpability, and

*AGY*recoding, indicating lumpability.

### Example 2 -- Primate nuclear DNA

*PP*-plot in Figure 7 clearly shows that the data are consistent with evolution under these conditions. Hence, it is appropriate to use our likelihood-ratio test to determine whether any of the recoding schemes would retain the Markovian properties of the original data. Figure 8 presents the

*PP*-plots from the likelihood-ratio test for lumpability for

*RY*-recoding, showing strong evidence against lumpability, and for the

*SW*-recoding, which provided the least evidence against lumpability. It is evident that no recoding should be applied to these data.

## Conclusions

Bias in estimates of phylogenetic parameters can occur when recoding of nucleotides or amino acids is used to transform data associated with models of evolution, which are not lumpable with respect to the recoding scheme used. A test proposed in this paper, which is based on a likelihood-ratio test, can yield an indication of whether the same results for estimable parameters can be expected from fitting a given model of evolution and its recoded version to the data.

## Declarations

### Acknowledgements

VAVR was supported by the University of Sydney World Scholars scheme. LSJ, KWL, and JR were supported by a Discovery Grant (DP0453173) from the Australian Research Council. We thank DL Lovell and TKF Wong for their constructive suggestions to this manuscript.

**Declarations**

The publication costs for this article were funded by resources made available to LSJ by CSIRO and JR by the School of Mathematics and Statistics, University of Sydney.

This article has been published as part of *BMC Bioinformatics* Volume 15 Supplement 2, 2014: Selected articles from the Twelfth Asia Pacific Bioinformatics Conference (APBC 2014): Bioinformatics. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/15/S2.

## Authors’ Affiliations

## References

- Irwin DM, Kocher TD, Wilson AC: Evolution of the cytochrome b gene in mammals. Journal of Molecular Evolution. 1991, 32: 128-144. 10.1007/BF02515385.View ArticlePubMedGoogle Scholar
- Adkins RM, Honeycutt RL: Molecular phylogeny of the superorder Arconta. Proceedings of the National Academy of Science of the United States of America. 1991, 88: 10317-10321. 10.1073/pnas.88.22.10317.View ArticleGoogle Scholar
- Adkins RM, Honeycutt RL: Evolution of the primate cytochrome c oxidase subunit II gene. Journal of Molecular Evolution. 1994, 38: 215-231.View ArticlePubMedGoogle Scholar
- Woese CR, Achenbach L, Rouviere P, Mandelco L: Archaeal phylogeny: reexamination of the phylogenetic position of
*Archaeoglobus fulgidus*in light of certain composition-induced artifacts. Systematic and Applied Microbiology. 1991, 14: 364-371. 10.1016/S0723-2020(11)80311-5.View ArticlePubMedGoogle Scholar - Phillips MJ, Penny D: The root of the mammalian tree inferred from whole mithocondrial genomes. Molecular Phylogenetics and Evolution. 2003, 28: 171-185. 10.1016/S1055-7903(03)00057-5.View ArticlePubMedGoogle Scholar
- Cavender JA, Felsenstein J: Invariants of phylogenies in a simple case with discrete states. Journal of Classification. 1987, 4: 57-71. 10.1007/BF01890075.View ArticleGoogle Scholar
- Gibson A, Gowri-Shankar V, Higgs PG, Rattray M: A comprehensive analysis of mammalian mithochondrial genome base composition and improved phylogenetic methods. Molecular Biology and Evolution. 2005, 22: 251-264.View ArticlePubMedGoogle Scholar
- Millen RS, Olmstead RG, Adams KL, Palmer JD, Lao NT, Heggie L, Kavanagh TA, Hibberd JM, Gray JC, Morden CW, Calie PJ, Jermiin LS, Wolfe KH: Many parallel losses of infA from chloroplast DNA during angiosperm evolution with multiple independent transfers to the nucleus. The Plant Cell. 2001, 13: 645-658. 10.1105/tpc.13.3.645.PubMed CentralView ArticlePubMedGoogle Scholar
- Phillips MJ, Lin YH, Harrison GL, Penny D: Mitochondrial genomes of a bandicoot and a brushtail possum confirm the monophyly of australidelphian marsupials. Proceedings of the Royal Society London Series B. 2001, 268: 1533-1538. 10.1098/rspb.2001.1677.View ArticleGoogle Scholar
- Kosiol C, Goldman N, Buttimore NH: A new criterion and method for amino acid classification. Journal of Theoretical Biology. 2004, 228: 97-106. 10.1016/j.jtbi.2003.12.010.View ArticlePubMedGoogle Scholar
- Kosiol C: Markov models for protein sequence evolution. PhD thesis. 2006, University of CambridgeGoogle Scholar
- Phillips MJ, Delsuc F, Penny D: Genome-scale phylogeny and the detection of systematic biases. Molecular Biology and Evolution. 2004, 21: 1455-1458. 10.1093/molbev/msh137.View ArticlePubMedGoogle Scholar
- Ho JWK, Adams CE, Lew JB, Matthews TJ, Ng CC, Shahabi-Sirjani A, Tan LH, Zhao Y, Easteal S, Wilson SR, Jermiin LS: SeqVis: Visualization of compositional heterogeneity in large alignments of nucleotides. Bioinformatics. 2006, 221: 2162-2163.View ArticleGoogle Scholar
- Susko E, Roger AJ: On reduced amino acid alphabets for phylogenetic inference. Molecular Biology and Evolution. 2007, 24: 2139-2150. 10.1093/molbev/msm144.View ArticlePubMedGoogle Scholar
- Anisimova M, Kosiol C: Investigating protein-coding sequence evolution with probabilistic codon substitution models. Molecular Biology and Evolution. 2004, 26: 255-271.View ArticleGoogle Scholar
- Masta SE, Longhorn SJ, Boore JL: Arachnid relationships based on mitochondrial genomes: asymmetric nucleotide and amino acid bias affects phylogenetic analyses. Molecular Phylogenetics and Evolution. 2009, 50: 117-128. 10.1016/j.ympev.2008.10.010.View ArticlePubMedGoogle Scholar
- Phillips MJ, Gibb GC, Crimp EA, Penny D: Tinamous and moa flock together: mitochondrial genome sequence analysis reveals independent losses of flight among ratites. Systematic Biology. 2010, 59: 90-107. 10.1093/sysbio/syp079.View ArticlePubMedGoogle Scholar
- Regier JC, Shultz JW, Zwick A, Hussey A, Ball B, Wetzer R, Martin JW, Cunningham CW: Arthropod relationships revealed by phylogenomic analysis of nuclear protein-coding sequences. Nature. 2010, 463: 1079-1083. 10.1038/nature08742.View ArticlePubMedGoogle Scholar
- Criscuolo A, Gribaldo S: Large-scale phylogenomic analyses indicate a deep origin of primary plastids within Cyanobacteria. Molecular Biology and Evolution. 2011, 28: 3019-3032. 10.1093/molbev/msr108.View ArticlePubMedGoogle Scholar
- Regier JC, Zwick A: Sources of signal in 62 protein-coding nuclear genes for higher-level phylogenetics of arthropods. PLoS ONE. 2011, 6: e23408-10.1371/journal.pone.0023408.PubMed CentralView ArticlePubMedGoogle Scholar
- Cho S, Zwick A, Regier JC, Mitter C, Cummings MP, Yao J, Du Z, Zhao H, Kawahara AY, Weller S, Davis DR, Baixeras J, Brown JW, Parr C: Can deliberately incomplete gene sample augmentation improve a phylogeny estimate for the advanced moths and butterflies (Hexapoda: Lepidoptera)?. Systematic Biology. 2011, 60: 782-796. 10.1093/sysbio/syr079.PubMed CentralView ArticlePubMedGoogle Scholar
- White NE, Phillips MJ, Gilbert MTP, Alfaro-Nunez A, Willerslev E, Mawson PR, Spencer PBS, Bunce M: The evolutionary history of cockatoos (Aves: Psittaciformes: Cacatuidae). Molecular Phylogenetics and Evolution. 2011, 59: 615-622. 10.1016/j.ympev.2011.03.011.View ArticlePubMedGoogle Scholar
- Zwick A, Regier JC, Cummings MP, Mitter C: Increased gene sampling yields robust support for higher-level clades within Bombycoidea (Lepidoptera). Systematic Entomology. 2011, 36: 31-43. 10.1111/j.1365-3113.2010.00543.x.View ArticleGoogle Scholar
- Niehuis O, Hartig G, Grath S, Pohl H, Lehmann J, Tafer H, Donath A, Krauss V, Eisenhardt C, Hertel J, Petersen M, Mayer C, Meusemann K, Peters RS, Stadler PF, Beutel RG, Bornberg-Bauer E, McKenna DD, Misof B: Genomic and morphological evidence converge to resolve the enigma of Strepsiptera. Current Biology. 2012, 22: 1309-1313. 10.1016/j.cub.2012.05.018.View ArticlePubMedGoogle Scholar
- Regier JC, Brown JW, Mitter C, Baixeras J, Cho S, Cummings MP, Zwick A: A molecular phylogeny for the leaf-roller moths (Lepidoptera: Tortricidae) and its implications for classification and life history evolution. PLoS ONE. 2012, 7: e35574-10.1371/journal.pone.0035574.PubMed CentralView ArticlePubMedGoogle Scholar
- Regier JC, Mitter C, Solis MA, Hayden JE, Landry B, Nuss M, Simonsen TJ, Yen S-H, Zwick A, Cummings MP: A molecular phylogeny for the pyraloid moths (Lepidoptera: Pyraloidea) and its implications for higher-level classification. Systematic Entomology. 2012, 37: 635-656. 10.1111/j.1365-3113.2012.00641.x.View ArticleGoogle Scholar
- Zwick A, Regier JC, Zwickl DJ: Resolving discrepancy between nucleotides and amino acids in deeplevel arthropod phylogenomics: differentiating serine codons in 21-amino-acid models. PLoS ONE. 2012, textbf7: e47450-View ArticleGoogle Scholar
- Gibb GC, Kennedy M, Penny D: Beyond phylogenetics and evolution: pelecaniform and Ciconiiform birds, and long-term niche stability. Molecular Phylogenetics and Evolution. 2013, 68: 229-238. 10.1016/j.ympev.2013.03.021.View ArticlePubMedGoogle Scholar
- Regier JC, Mitter C, Zwick A, Bazinet AL, Cummings MP, Kawahara AY, Sohn J-C, Zwickl DJ, Cho S, Davis DR, Baixeras J, Brown J, Parr C, Weller S, Lees DC, Mitter KT: A large-scale, higher-level, molecular phylogenetic study of the insect order Lepidoptera (Moths and Butterflies). PLoS ONE. 2013, 8: e58568-10.1371/journal.pone.0058568.PubMed CentralView ArticlePubMedGoogle Scholar
- Rota-Stabelli O, Lartillot N, Philippe H, Pisani D: Serine codon-usage bias in deep phylogenomics: pancrustacean relationships as a case study. Systematic Biology. 2013, 62: 121-133. 10.1093/sysbio/sys077.View ArticlePubMedGoogle Scholar
- Sohn J-C, Regier JC, Mitter C, Davis D, Landry J-F, Zwick A, Cummings MP: A molecular phylogeny for Yponomeutoidea (Insecta, Lepidoptera, Ditrysia) and its implications for classification, biogeography and the evolution of host plant use. PLoS ONE. 2013, textbf8: e55066-View ArticleGoogle Scholar
- Lau KW: Studies of methods used to infer molecular phylogeny: Dealing with the effect of compositional heterogeneity. PhD thesis. 2009, University of Sydney, School of Biological Sciences;Google Scholar
- Guédon Y, d'Aubenton-Carafa Y, Thermes C: Analysing grouping of nucleotides in DNA sequences using lumped processes constructed from Markov chains. Journal of Mathematical Biology. 2006, 52: 343-372. 10.1007/s00285-005-0358-y.View ArticlePubMedGoogle Scholar
- Swofford DL, Olsen GJ, Waddell PJ, Hillis DM: Phylogenetic inference. Molecular Systematics. Edited by: Hillis DM, Moritz C, Mable BK. 1996, Sunderland: Sinauer Associates, 407-514.Google Scholar
- Nomenclature Committee of the International Union of Biochemistry, (NC-IUB): Nomenclature for Incompletely Specified Bases in Nucleic Acid Sequences: Recommendations 1984. Proceedings of the National Academy of Sciences of the United States of America. 1986, 83: 4-8.View ArticleGoogle Scholar
- Bryant D, Galtier N, Poursat MA: Likelihood calculation in molecular phylogenetics. Mathematics evolution and phylogeny. Edited by: Gascuel O. 2005, New York: Oxford University Press, 33-92.Google Scholar
- Jayaswal V, Jermiin LS, Robinson J: Estimation of phylogeny using a general Markov model. Evolutionary Bioinformatics. 2005, 1: 62-80.Google Scholar
- Ababneh F, Jermiin LS, Ma C, Robinson J: Matched-pairs tests of homogeneity with applications to homologous nucleotide sequences. Bioinformatics. 2006, 22: 1225-1231. 10.1093/bioinformatics/btl064.View ArticlePubMedGoogle Scholar
- Jermiin LS, Jayaswal V, Ababneh F, Robinson J: Phylogenetic model evaluation. Bioinformatics: Data, sequence analysis, and evolution − Volume 1. Edited by: Keith J. 2008, Humana Press. Totawa, 331-363.View ArticleGoogle Scholar
- Ababneh F, Jermiin LS, Robinson J: Generation of the exact distribution and simulation of matched nucleotide sequences on a phylogenetic tree. Journal of Mathematical Modelling and Algorithms. 2006, 5: 291-303. 10.1007/s10852-005-9017-y.View ArticleGoogle Scholar
- Iosifescu M: Finite Markov processes and their applications. 1980, Chichester: John Wiley and Sons, LtdGoogle Scholar
- Kemeny JG, Snell JL: Finite Markov chains. 1983, New York: Springer-VerlagGoogle Scholar
- Jukes TH, Cantor CR: Evolution of protein molecules. Mammalian Protein Metabolism. Edited by: Munro HN. 1969, Academic Press. New York, 21-132.View ArticleGoogle Scholar
- Felsenstein J: Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of Molecular Evolution. 1981, 17: 368-376. 10.1007/BF01734359.View ArticlePubMedGoogle Scholar
- Kimura M: A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. Journal of Molecular Evolution. 1980, 16: 111-120. 10.1007/BF01731581.View ArticlePubMedGoogle Scholar
- Kishino H, Yano T: Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution. 1985, 22: 160-174. 10.1007/BF02101694.View ArticlePubMedGoogle Scholar
- Jernigan RW, Baran RH: Testing lumpability in Markov chains. Statistics and Probability Letters. 2003, 64: 17-23. 10.1016/S0167-7152(03)00126-3.View ArticleGoogle Scholar
- Schweder T, Spjotvoll E: Plots of P-values to evaluate many tests simultaneously. Biometrika. 1982, 69: 493-502. 10.1093/biomet/69.3.493.View ArticleGoogle Scholar
- Lanave C, Preparata G, Saccone C, Serio G: A new method for calculating evolutionary substitution rates. Journal of Molecular Evolution. 1984, 20: 86-93. 10.1007/BF02101990.View ArticlePubMedGoogle Scholar
- Tosi AJ, Detwiler KM, Disotell TR: X-chromosomal window into the evolutionary history of the guenons (Primates: Cercopithecini). Molecular Phylogenetics and Evolution. 2005, 36: 58-66. 10.1016/j.ympev.2005.01.009.View ArticlePubMedGoogle Scholar
- Benson DA, Cavanaugh M, Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, Sayers EW: Genbank. Nucleic Acids Research. 2013, 41: D36-D42. 10.1093/nar/gks1195.PubMed CentralView ArticlePubMedGoogle Scholar
- Katoh K, Standley DM: MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Molecular Biology and Evolution. 2013, 30: 772-780. 10.1093/molbev/mst010.PubMed CentralView ArticlePubMedGoogle Scholar
- Gouy M, Guindon S, Gascuel O: SeaView version 4: a multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Molecular Biology and Evolution. 2010, 27: 221-224. 10.1093/molbev/msp259.View ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.