Impossibility of consistent distance estimation from sequence lengths under the TKF91 model

Phylogeny estimation consists in the inference of an evolutionary tree from extant species data, commonly molecular sequences (e.g. DNA, amino acid). A large body of theoretical work exists on the statistical properties of standard reconstruction methods [Ste16, War17]. Typically in such analyses, one assumes that sequences have evolved on a fixed rooted tree, from a common ancestor sequence to the leaf sequences, according to some Markovian stochastic process. Often these processes model site substitutions exclusively, with the underlying assumption being that the data has been properly aligned in a pre-processing step. In contrast, relatively little theoretical work has focused on models of insertions and deletions (indels) together with substitutions, in spite of the fact that such models have been around for some time [TKF91, TKF92]. See e.g. [Tha06, DR13, ARS15, FR20].

One extra piece of information available under indel models is the length of the sequence, which itself evolves according to a Markov process on the tree. The notable work of Thatte [Tha06] shows that leaf sequence lengths alone are in fact enough to reconstruct phylogenies, through a distance-based approach. More specifically, it is shown in [Tha06, (27)] that under the TKF91 model [TKF91] the expectation of the sequence length $N_{v}$ at a leaf $v$ conditioned on the sequence length $N_{u}$ at another leaf $u$ separated from $v$ by an amount of time $t_{uv}$ is

$\displaystyle\mathcal{N}_{v}(t)=\bar{L}+\left(N_{u}-\bar{L}\right)e^{-\mu t_{uv}(1-\lambda/\mu)}$

(1)

where $\bar{L}=\frac{\lambda/\mu}{1-\lambda/\mu}$ is the expected length at stationarity, where $\lambda<\mu$ are the rates of insertion and deletion respectively. (Full details on the TKF91 model are provided in Section 2.) Hence we see from (1) that the full distribution of sequence lengths suffices to recover $\lambda/\mu$ and all $\mu t_{uv}$’s.

The tree topology can then be recovered using standard results about the metric properties of phylogenies [Ste16]. That is, the tree is identifiable from the sequence lengths under the TKF91 model in the sense that two distinct tree topologies $T_{1}\neq T_{2}$ necessarily produce distinct joint distributions of sequence lengths at the leaves.

It is also suggested in [Tha06]—without a full rigorous proof—that the scheme above could be used to reconstruct phylogenies from a single sample of sequence lengths at the leaves in the limit where $\lambda\nearrow\mu$. The latter asymptotics ensure that the expected sequence length at stationarity $\bar{L}$ diverges, and serves as a proxy for the amount of data growing to infinity. However, in this short note, we show that no consistent distance estimator exists in this limit. Formally we establish that the distributions of pairs of sequence lengths at different distances cannot be distinguished with probability going to $1$ as $\lambda\nearrow\mu$. Hence, while the tree is identifiable from the distribution of the sequence lengths at the leaves, one sample of the sequence lengths alone cannot be used in a distance-based approach of the type described above to reconstruct the tree consistently as $\lambda\nearrow\mu$. On the technical side our proof follows by noting that, under the TKF91 model, the sequence length is (morally) a sum of independent random variables with finite variances, to which we apply a central limit theorem. One complication is to obtain a limit theorem that is uniform in the parameter $\lambda/\mu$. We expect that our techniques will be useful to analyze other bioinformatics methods under indel processes, for instance methods based on $k$-mer statistics (see e.g. [YZ08, Hau13]). Further intuition on our results is provided in Section 3.

Organization. The rest of the paper is organized as follows. The TKF91 model is reviewed in Section 2. Our main result, together with a proof sketch, is stated in Section 3. Details of the proof are provided in Section 4.

In this section, we recall the TKF91 sequence evolution model [TKF91]. To simplify the presentation, we restrict ourselves to a two-state version of the model, as we will only require the underlying sequence-length process.

We will be concerned with the underlying sequence-length process.

We are interested in this process on a rooted tree $T$. Denote the index set by $\Gamma_{T}$. The root vertex $\rho$ is assigned a state $\mathcal{I}_{\rho}\in\mathcal{S}$, drawn from stationary distribution on $\mathcal{S}$. This state is then evolved down the tree according to the following recursive process. Moving away from the root, along each edge $e=(u,v)\in E$, conditionally on the state $\mathcal{I}_{u}$, we run the indel process for a time $\ell_{(u,v)}$. Denote by $\mathcal{I}_{t}$ the resulting state at $t\in e$. Then the full process is denoted by $\{\mathcal{I}_{t}\}_{t\in\Gamma_{T}}$. In particular, the set of leaf states is $\mathcal{I}_{\partial T}=\{\mathcal{I}_{v}:v\in\partial T\}$.

Setting. Throughout this paper, we let $\mathbb{P}_{\vec{x}}$ be the probability measure when the root state is $\vec{x}$. If the root state is chosen according to a distribution $\nu$, then we denote the probability measure by $\mathbb{P}_{\nu}$. We also denote by $\mathbb{P}_{M}$ the conditional probability measure for the event that the root state has length $M$.

For our purposes, it will suffice to focus on the space $\mathcal{T}_{2}$ of star trees with two leaves that have the same finite distance $h\in(0,\infty)$ from the root and are labeled as $\{1,2\}$. This distance $h$ is the height of the tree. The indel process on a tree $T\in\mathcal{T}_{2}$ reduces to a pair of indel processes $(\mathcal{I}_{t}^{1},\mathcal{I}_{t}^{2})_{t\geq 0}$ that are independent upon conditioning on the root state $\mathcal{I}_{\rho}=\mathcal{I}_{0}^{1}=\mathcal{I}_{0}^{2}$. We always assume the root state is chosen according to the equilibrium distribution $\Pi$. So the distribution of $(\mathcal{I}_{0}^{1},\mathcal{I}_{0}^{2})\in\mathcal{S}\times\mathcal{S}$ is

$$\widehat{\nu}_{0}(\vec{x},\vec{y})=\begin{cases}\Pi(\vec{x})&\text{ if }\vec{x}=\vec{y},\\ 0&\textnormal{otherwise.}\end{cases}$$

Our main theorem is an impossibility result: the distributions of pairs of sequence lengths at different distances cannot be distinguished with probability going to $1$ as $\lambda\nearrow\mu$. Following [Tha06], we consider the asymptotic regime where $\lambda\nearrow\mu$, which implies that the expected sequence length at stationarity tends to $+\infty$. Recall that the total variation distance between two probability measures $\tau_{1}$ and $\tau_{2}$ on a countable measure space $E$ is

$\displaystyle\|\tau_{1}-\tau_{2}\|_{TV}=\frac{1}{2}\sum_{\sigma\in E}\left|\tau_{1}(\sigma)-\tau_{2}(\sigma)\right|.$

Recall that the total variation distance can be written as

$\displaystyle\|\tau_{1}-\tau_{2}\|_{TV}=\sup_{A\subseteq E}\left|\tau_{1}(A)-\tau_{2}(A)\right|.$

So (3) implies that there is no test that can distinguish between $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ with probability going to $1$ as $\lambda\nearrow\mu$.

Proof idea. We first give a heuristic argument that underlies our formal proof. Without loss of generality, assume that the deletion rate is $\mu=1$. The stationary length $M$ at the root is geometric with mean and standard deviation both of order $1/(1-\lambda)$. So we can think of the root length as roughly $M\approx C/(1-\lambda)$ with significant probability. Ignoring the small effect of the immortal link and conditioning on $M$, the lengths at the leaves are sums of independent random variables, specifically the progenies of the $M$ mortal links of the root. The mean and variance of these variables can be computed explicitly from continuous-time Markov chain theory (see (10) below; see also [Tha06, (27), (31)]). As $\lambda\nearrow 1$, the difference in expectation between heights $h_{1}$ and $h_{2}$ turns out to be

$\displaystyle Me^{-(1-\lambda)h_{1}}-Me^{-(1-\lambda)h_{2}}\approx\frac{C}{1-\lambda}[-(1-\lambda)h_{1}+(1-\lambda)h_{2}]\approx C(h_{2}-h_{1}),$

(4)

while the variance is of order

$\displaystyle M\frac{e^{-(1-\lambda)h_{i}}(1-e^{-(1-\lambda)h_{i}})}{1-\lambda}\approx\frac{C}{1-\lambda}\frac{(1-\lambda)h_{i}}{1-\lambda}\approx\frac{Ch_{i}}{1-\lambda}.$

(5)

The key observation is that the variance (5) is $\gg$ than the square of the expectation difference (4). Hence, by the central limit theorem, one can expect significant overlap between the length distributions under $h_{1}$ and $h_{2}$, making them hard to distinguish even as $\lambda\nearrow 1$. We formalize this argument next.

We observe that (3) is equivalent to

$$\liminf_{\lambda\nearrow 1}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\mathbb{P}_{\Pi}(\vec{N}^{(1)}=\vec{y})\wedge\mathbb{P}_{\Pi}(\vec{N}^{(2)}=\vec{y})>0.$$

(6)

Indeed the total variation distance between two probability measures $\tau_{1}$ and $\tau_{2}$ on a countable space $E$ can also be written as

$\displaystyle\|\tau_{1}-\tau_{2}\|_{TV}=1-\sum_{\sigma\in E}\tau_{1}(\sigma)\land\tau_{2}(\sigma).$

The rest of the proof is to establish (6). It involves a series of steps:

1. 1.

   We first reduce the problem to a sum of independent random variables by conditioning on the root sequence length and ignoring the immortal link. In particular, we use the fact that there is a fairly uniform probability that $M$ is in an interval of size $1/(1-\lambda)$ around $1$. And we remove the effect of the immortal link by conditioning on its having no descendant, an event of positive probability.

2. 2.

   The central limit theorem (CLT) implies that there is a significant overlap between the two sums. More precisely, we need a local CLT (see e.g. [Dur10]) to derive the sort of pointwise lower bound needed in (6). However the bound we require must be uniform in $\lambda$ and we did not find in the literature a result of quite this form. Instead, we use an argument based on the Berry-Esséen theorem (again see e.g. [Dur10]). We first establish overlap over $\Omega(\sqrt{M})$ constant size intervals for the sum of the first $M-1$ mortal links, and then we use the final mortal link to match the probabilities on common point values under heights $h_{1}$ and $h_{2}$.

3. 3.

   Finally we bound the sum in (6).

In this section we give the details of the proof of Theorem 1. We follow the steps described in the previous section.

Precisely, for any $\lambda_{*}\in(0,1)$ and $0<c_{1}<1<c_{2}<+\infty$, using (2)

	$\displaystyle\liminf_{\lambda\nearrow 1}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\mathbb{P}_{\Pi}(\vec{N}^{(1)}=\vec{y})\wedge\mathbb{P}_{\Pi}(\vec{N}^{(2)}=\vec{y})$
	$\displaystyle\geq\inf_{\lambda\in(\lambda_{*},1)}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\mathbb{P}_{\Pi}(\vec{N}^{(1)}=\vec{y})\wedge\mathbb{P}_{\Pi}(\vec{N}^{(2)}=\vec{y})$
	$\displaystyle=\inf_{\lambda\in(\lambda_{*},1)}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\left[\sum_{M\in\mathbb{Z}_{+}}\gamma^{(\lambda)}_{M}\,\mathbb{P}_{M}(\vec{N}^{(1)}=\vec{y})\right]\wedge\left[\sum_{M\in\mathbb{Z}_{+}}\gamma^{(\lambda)}_{M}\,\mathbb{P}_{M}(\vec{N}^{(2)}=\vec{y})\right]$
	$\displaystyle=\inf_{\lambda\in(\lambda_{*},1)}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\sum_{M\in\mathbb{Z}_{+}}(1-\lambda)\lambda^{M}\left[\mathbb{P}_{M}(\vec{N}^{(1)}=\vec{y})\wedge\mathbb{P}_{M}(\vec{N}^{(2)}=\vec{y})\right]$
	$\displaystyle\geq\inf_{\lambda\in(\lambda_{*},1)}\sum_{M\in\left[\frac{c_{1}}{1-\lambda},\frac{c_{2}}{1-\lambda}\right]}(1-\lambda)\lambda^{M}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\mathbb{P}_{M}(\vec{N}^{(1)}=\vec{y})\wedge\mathbb{P}_{M}(\vec{N}^{(2)}=\vec{y})$
	$\displaystyle\geq c_{3}(c_{2}-c_{1})\inf_{\lambda\in(\lambda_{*},1)}\inf_{M\in\left[\frac{c_{1}}{1-\lambda},\frac{c_{2}}{1-\lambda}\right]}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\mathbb{P}_{M}(\vec{N}^{(1)}=\vec{y})\wedge\mathbb{P}_{M}(\vec{N}^{(2)}=\vec{y}),$		(7)

where $c_{3}$ is a lower bound on $\lambda^{M}$ for $M\in\left[\frac{c_{1}}{1-\lambda},\frac{c_{2}}{1-\lambda}\right]$ and $\lambda\in(\lambda_{*},1)$.

Let $\mathcal{Z}_{0}$ be the event that the immortal link of the root sequence produces no mortal link in either leaf sequences. Let $\mathbb{P}_{M,\bullet}$ be the probability conditioned on that event, and $c_{4}$ be a lower bound on the probability of $\mathcal{Z}_{0}$ uniform in $\lambda\in(\lambda_{*},1)$. Under $\mathbb{P}_{M,\bullet}$, the two components of $\vec{N}^{(1)}$ are conditionally independent and each is a sum of $M$ i.i.d. random variables corresponding to the progenies of mortal links. Hence (7) is at least

	$\displaystyle c_{4}c_{3}(c_{2}-c_{1})\inf_{\lambda\in(\lambda_{*},1)}\inf_{M\in\left[\frac{c_{1}}{1-\lambda},\frac{c_{2}}{1-\lambda}\right]}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\mathbb{P}_{M,\bullet}(\vec{N}^{(1)}=\vec{y})\wedge\mathbb{P}_{M,\bullet}(\vec{N}^{(2)}=\vec{y})$
	$\displaystyle\geq c_{4}c_{3}(c_{2}-c_{1})\inf_{\lambda\in(\lambda_{*},1)}\inf_{M\in\left[\frac{c_{1}}{1-\lambda},\frac{c_{2}}{1-\lambda}\right]}\sum_{\vec{y}\in\mathbb{Z}_{+}^{2}}\left[p_{M,y_{1}}^{(\lambda)}(h_{1})\,p_{M,y_{2}}^{(\lambda)}(h_{1})\right]\wedge\left[p_{M,y_{1}}^{(\lambda)}(h_{2})\,p_{M,y_{2}}^{(\lambda)}(h_{2})\right].$		(8)

where we let $p_{y_{1},y_{2}}^{(\lambda)}(t)=\mathbb{P}_{y_{1},\bullet}(|\mathcal{I}_{t}|=y_{2})$ be the transition probability of the length process without the immortal link.

The sum in (8) leads us to study the overlap between the probability distributions $p_{M,\,\cdot}^{(\lambda)}(t):=\{p_{M,j}^{(\lambda)}(t)\}_{j\in\mathbb{Z}_{+}}$ for $t=h_{1},h_{2}$ and $M\in\left[\frac{c_{1}}{1-\lambda},\frac{c_{2}}{1-\lambda}\right]$. The central limit theorem is what we need. However, because of our need for a bound that is uniform in $\lambda$, we shall apply the Berry-Esséen theorem. Specifically, we apply the latter bound to the progenies of the first $M-1$ mortal links of the root sequence. The idea is to show that $\Omega(\sqrt{M})$ summands in (8) have value $\Omega(1/\sqrt{M})$, for each of $h_{1}$ and $h_{2}$ separately, and then use the last mortal link to “match” all these values between $h_{1}$ and $h_{2}$.

Let the mean and the variance of $L^{i}_{t}$ be

$$\beta:=\beta(\lambda,t):={\mathbb{E}}[L^{i}_{t}]\quad\text{and}\quad\sigma^{2}:=\sigma^{2}(\lambda,t):={\mathbb{E}}|L^{i}_{t}-\beta|^{2}.$$

(9)

As is expected, the distribution $p_{M,\,\cdot}^{(\lambda)}(t)$ is approximately Gaussian with mean $\beta M$ and variance $\sigma^{2}M$. We quantify this statement in the bound (11) below, after proving some moment bounds.

Let $F_{M}^{(\lambda)}(t)$ be the cumulative distribution function (CDF) of the probability distribution $p_{M,\cdot}^{(\lambda)}(t)$. That is,

$$F_{M}^{(\lambda)}(t)(x)=\sum_{j:j\leq x}p_{M,j}^{(\lambda)}(t)={\mathbb{P}}(S_{M}(t)\leq x).$$

Now consider the interval with length roughly the standard deviation and centered at around one of the means, $\beta(\lambda,h_{1})(M-1)$. Then consider an equi-partition of this interval into roughly $\sqrt{M-1}$ many pieces of constant length. Precisely, we write $\beta_{1}:=\beta(\lambda,h_{1})$ and $\sigma_{1}:=\sigma(\lambda,h_{1})$ for simplicity. Then for an arbitrary constant $K>0$,

$$\Big{(}\beta_{1}(M-1)-\sigma_{1}\sqrt{M-1},\;\beta_{1}(M-1)+\sigma_{1}\sqrt{M-1}\Big{)}=\bigcup_{r\in\Lambda^{M}_{K}}\mathcal{J}^{M}_{r}(c)$$

where the sub-intervals

$$\mathcal{J}^{M}_{r}(K):=\big{(}\beta_{1}(M-1)+r\sigma_{1}K,\;\beta_{1}(M-1)+(r+1)\sigma_{1}K\big{)}$$

have constant width $\sigma_{1}K$ for $r\in\Lambda^{M}_{K}$, and

$$\Lambda^{M}_{K}:=\left\{-\sigma_{1}\sqrt{M-1},\,-\sigma_{1}(\sqrt{M-1}-K)\,\ldots,0,\,\sigma_{1}K,\,\ldots,\sigma_{1}(\sqrt{M-1}-K)\right\}.$$

Lemma 3 below says that there exists a constants $c=c_{7}$ and $K$ large enough (depending on $c_{5}$ and $c_{6}$ but not on $\lambda$) such that each of these intervals contains mass at least $\frac{c_{8}}{\sqrt{M-1}}$ under both probability distributions $p_{M-1,\,\cdot}^{(\lambda)}(h_{1})$ and $p_{M-1,\,\cdot}^{(\lambda)}(h_{2})$. See Figure 1. Write $p_{M-1,\,A}^{(\lambda)}(t)=\sum_{j\in A}p_{M-1,\,j}^{(\lambda)}(t)$ for simplicity.

SR was supported by NSF grants DMS-1614242, CCF-1740707 (TRIPODS), DMS-1902892, and DMS-1916378, as well as a Simons Fellowship and a Vilas Associates Award. BL was supported by DMS-1614242, CCF-1740707 (TRIPODS), DMS-1902892 (to SR). WTF was supported by NSF grant DMS-1614242 (to SR) and DMS-1855417.