Julia A Palacios
email: juliapr@stanford.edu

Amandine Véber

Lorenzo Cappello

Zhangyuan Wang

John Wakeley

Sohini Ramachandran

Bayesian Estimation of Population Size Changes by Sampling Tajima's Trees

The large state space of gene genealogies is a major hurdle for inference methods based on Kingman's coalescent. Here, we present a new Bayesian approach for inferring past population sizes which relies on a lower resolution coalescent process we refer to as "Tajima's coalescent".

Tajima's coalescent has a drastically smaller state space, and hence it is a computationally more efficient model, than the standard Kingman coalescent. We provide a new algorithm for efficient and exact likelihood calculations for data without recombination, which exploits a directed acyclic graph and a correspondingly tailored Markov Chain Monte Carlo method. We compare the performance of our Bayesian Estimation of population size changes by Sampling Tajima's Trees (BESTT) with a popular implementation of coalescent-based inference in BEAST using simulated data and human data. We empirically demonstrate that BESTT can accurately infer effective population sizes, and it further provides an efficient alternative to the Kingman's coalescent. The algorithms described here are implemented in the R package phylodyn, which is available for download at https://github.com/JuliaPalacios/phylodyn.

Introduction

Modeling gene genealogies from an alignment of sequences -timed and rooted bifurcating trees reflecting the ancestral relationships among sampled sequences -is a key step in coalescent-based inference of evolutionary parameters such as effective population sizes. In the neutral coalescent model without recombination, observed sequence variation is produced by a stochastic process of mutation acting along the branches of the gene genealogy (Kingman, 1982a;[START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF], which is modeled as a realization of the coalescent point process at a neutral non-recombining locus. In the coalescent point process, the rate of coalescence (the merging of two lineages into a common ancestor at some time in the past) is a function that varies with time, and it is inversely proportional to the effective population size at time t, N (t) (Kingman, 1982b;[START_REF] Slatkin | Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations[END_REF] without recombination [START_REF] Griffiths | Monte Carlo inference methods in population genetics[END_REF][START_REF] Kuhner | Comparing likelihood and Bayesian coalescent estimation of population parameters[END_REF][START_REF] Minin | Smooth skyride through a rough skyline: Bayesian coalescent-based inference of population dynamics[END_REF][START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF][START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF][START_REF] Palacios | Gaussian process-based Bayesian nonparametric inference of population trajectories from gene genealogies[END_REF][START_REF] Gill | Improving bayesian population dynamics inference: A coalescent-based model for multiple loci[END_REF]. Our approach is of this type. Still other methods account for recombination across larger genomic segments [START_REF] Li | Inference of human population history from individual whole-genome sequences[END_REF][START_REF] Sheehan | Estimating variable effective population sizes from multiple genomes: A sequentially Markov conditional sampling distribution approach[END_REF][START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF][START_REF] Palacios | Bayesian nonparametric inference of population size changes from sequential genealogies[END_REF]. In spite of their variety, all these methods must contend with two major challenges: (i) choosing a prior distribution or functional form for (N (t)) t≥0 , and (ii) integrating over the large hidden state space of genealogies. For example, several previous approaches have assumed exponential growth [START_REF] Griffiths | Monte Carlo inference methods in population genetics[END_REF][START_REF] Kuhner | Maximum likelihood estimation of population growth rates based on the coalescent[END_REF][START_REF] Kuhner | Comparing likelihood and Bayesian coalescent estimation of population parameters[END_REF], in which case the estimation of (N (t)) t≥0 is reduced to the estimation of one or two parameters. In general, the functional form of (N (t)) t≥0 is unknown and needs to be inferred. A commonly used naive nonparametric prior on (N (t)) t≥0 is a piecewise linear or constant function defined on time intervals of constant or varying sizes [START_REF] Heled | Bayesian inference of population size history from multiple loci[END_REF][START_REF] Sheehan | Estimating variable effective population sizes from multiple genomes: A sequentially Markov conditional sampling distribution approach[END_REF][START_REF] Schiffels | Inferring human population size and separation history from multiple genome sequences[END_REF].

The specification of change points in such time-discretized effective population size trajectories is inherently difficult because it can lead to runaway behavior or large uncertainty in ( N (t)) t≥0 .

These difficulties can be avoided by the use of Gaussian-process priors in a Bayesian nonparametric framework, allowing accurate and precise estimation [START_REF] Palacios | Gaussian process-based Bayesian nonparametric inference of population trajectories from gene genealogies[END_REF][START_REF] Gill | Improving bayesian population dynamics inference: A coalescent-based model for multiple loci[END_REF][START_REF] Lan | An efficient Bayesian inference framework for coalescent-based nonparametric phylodynamics[END_REF][START_REF] Palacios | Bayesian nonparametric inference of population size changes from sequential genealogies[END_REF]. More precisely, the autocorrelation modeled with the Gaussian process avoids runaway behavior and large uncertainty in ( N (t)) t≥0 .

The second challenge for coalescent-based inference of (N (t)) t≥0 is the integration over the hidden state space of genealogies for large sample sizes. Given molecular sequence data Y at a single non-recombining locus and a mutation model with vector of parameters µ, current methods rely on calculating the marginal likelihood function Pr(Y|(N (t)) t≥0 , µ) by integrating over all possible coalescence and mutation events. Under the infinite-sites mutation model without intra-locus recombination [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF], this integration requires a computationally expensive importance sampling technique or Markov Chain Monte Carlo (MCMC) techniques (Griffiths and Tavaré, 1994a;[START_REF] Stephens | Inference in molecular population genetics[END_REF][START_REF] Hobolth | Importance sampling for the infinite sites model[END_REF][START_REF] Wu | Exact computation of coalescent likelihood for panmictic and subdivided populations under the infinite sites model[END_REF][START_REF] Gronau | Bayesian inference of ancient human demography from individual genome sequences[END_REF]. Moreover, a maximum likelihood estimate of (N (t)) t≥0 cannot be explicitly obtained; instead, it is obtained by exploring a grid of parameter values [START_REF] Tavaré | Part I: Ancestral inference in population genetics[END_REF]. For finite-sites mutation models, current methods approximate the marginal likelihood function by integrating over all possible genealogies via MCMC methods (Equation [START_REF]Z 8 -Z 10 -Z 3 -Z 6 -Z 4 -Z 7 -Z 5 -Z 4 -Z 6 -Z 1 -Z 2 -Z 1 -Z 2 -Z 0 2. Z 8 -Z[END_REF]; [START_REF] Kuhner | LAMARC 2.0: maximum likelihood and bayesian estimation of population parameters[END_REF]; [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF]). In both cases, the marginal likelihood may be expressed as Pr Y|(N (t)) t≥0 , µ = Pr(Y|g, µ)Pr g|(N (t)) t≥0 dg, [START_REF]Z 8 -Z 10 -Z 3 -Z 6 -Z 4 -Z 7 -Z 5 -Z 4 -Z 6 -Z 1 -Z 2 -Z 1 -Z 2 -Z 0 2. Z 8 -Z[END_REF] in which Pr(•) is used to denote both the probability of discrete variables and the density of continuous variables. The integral above involves an (n -1)-dimensional integral over n -1 coalescent times and a sum over all possible tree topologies with n leaves. Therefore, these methods require a very large number of MCMC samples, and exploration of the posterior space of genealogies con-tinues to be an active area of research [START_REF] Kuhner | Maximum likelihood estimation of population growth rates based on the coalescent[END_REF][START_REF] Rannala | Bayes estimation of species divergence times and ancestral population sizes using dna sequences from multiple loci[END_REF][START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF][START_REF] Whidden | Quantifying mcmc exploration of phylogenetic tree space[END_REF][START_REF] Aberer | An efficient independence sampler for updating branches in bayesian markov chain monte carlo sampling of phylogenetic trees[END_REF].

Current methods rely on the Kingman n-coalescent process to model the sample's ancestry.

However, the state space of genealogical trees grows superexponentially with the number of samples, making inference computationally challenging for large sample sizes. In this study, we develop a Bayesian nonparametric model that relies on Tajima's coalescent, a lower resolution coalescent process with a drastically smaller state space than that of Kingman's coalescent. In particular, we approximate the posterior distribution Pr((N (t)) t≥0 , g T , τ | Y, µ), where g T corresponds to the Tajima's genealogy of the sample (see Figure 1A and Section 2.4), (log N (t)) t≥0 has Gaussian process prior with precision hyperparameter τ that controls the degree of regularity, and mutations occur according to the infinite-sites model of [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF]. This results in a new efficient method for inferring (N (t)) t≥0 called Bayesian Estimation by Sampling Tajima's Trees (BESTT), with a drastic reduction in the state space of genealogies. We show using simulated data that BESTT can accurately infer effective population size trajectories and that it provides a more efficient alternative than Kingman's coalescent models.

Next, we start with an overview of BESTT, detail our representation of molecular sequence data and define the Tajima coalescent process. We then introduce a new augmented representation of sequence data as a directed acyclic graph (DAG). This representation allows us to both calculate the conditional likelihood under the Tajima coalescent model, and to sample tree topologies compatible with the observed data. We then provide an algorithm for likelihood calculations and develop an MCMC approach to efficiently explore the state space of unknown parameters. Finally, we compare our method to other methods implemented in BEAST [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF] and estimate the effective population size trajectory from human mtDNA data. We close with a discussion of possible extensions and limitations of the proposed model and implementation.

2 Methods/Theory

Overview of BESTT

Our objective in the implementation of BESTT is to estimate the posterior distribution of model parameters by replacing Kingman's genealogy with Tajima's genealogy g T . A Tajima's genealogy does not include labels at the tips (Figure 1): we do not order individuals in the sample but label only the lineages that are ancestral to at least two individuals (that is, we only label the internal nodes of the genealogy). Replacing Kingman's genealogy by Tajima's genealogy in our posterior distribution exponentially reduces the size of the state space of genealogies (Figure 1B). In order to compute Pr(Y|g T , µ), the conditional likelihood of the data conditioned on a Tajima's genealogy, we assume the infinite sites model of mutations and leverage a directed acyclic graph (DAG) representation of sequence data and genealogical information. Note that the overall likelihood, Eq. ( 1), will differ only by a combinatorial factor from the corresponding likelihood under the Kingman coalescent. Our DAG represents the data with a gene tree (Griffiths and Tavaré, 1994a), constructed via a modified version of the perfect phylogeny algorithm of [START_REF] Gusfield | Efficient algorithms for inferring evolutionary trees[END_REF]. This provides an economical representation of the uncertainty and conditional independences induced by the model and the observed data.

Under the infinite-sites mutation model, there is a one-to-one correspondence between observed sequence data and the gene tree of the data [START_REF] Gusfield | Efficient algorithms for inferring evolutionary trees[END_REF] the gene tree representation with the allocation of the number of observed mutations along the Tajima's genealogy to generate a DAG (Section 2.5). The conditional likelihood Pr(Y | g T , µ) is then calculated via a recursive algorithm that exploits the auxiliary variables defined in the DAG nodes, marginalizing over all possible mutation allocations (Section 2.6). We approximate the joint posterior distribution Pr((N (t)) t≥0 , g T , τ | Y, µ) via an MCMC algorithm using Hamiltonian Monte Carlo for sampling the continuous parameters of the model and a novel Metropolis-Hastings algorithm for sampling the discrete tree space.

Summarizing sequence data Y as haplotypes and mutation groups

Let the data consist of n fully linked haploid sequences or alignments of nucleotides at s segregating sites sampled from n individuals at time t = 0 (the present). Note that any labels we afix to the individuals are arbitrary in the sense that they will not enter into the calculation of the likelihood.

We further assume the infinite sites mutation model of [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF] with mutation parameter µ and known ancestral states for each of the sites. Then we can encode the data into a binary matrix Y of n rows and s columns with elements y i,j ∈ {0, 1}, where 0 indicates the ancestral allele.

In order to calculate the Tajima's conditional likelihood Pr(Y | g T , µ), we first record each haplotype's frequency and group repeated columns to form mutation groups; a mutation group corresponds to a shared set of mutations in a subset of the sampled individuals. We record the cardinality of each mutation group (i.e., the number of columns that show each mutation group).

In Figure 2A, there are two columns labeled "b", corresponding to two segregating sites which have the exact same pattern of allelic states across the sample. Further, two individuals carry the derived allele of mutation group "b", so in this case the frequency of haplotype 7 and the cardinality of mutation group "b" are both equal to 2. Likewise, haplotype 4 has frequency 1 and carries 5 mutations that are split into mutation groups "a", "f" and "g" (the latter is not shown in Figure 2A, but appears in Figure 2B) of respective cardinalities 1, 3 and 1. We denote the number of haplotypes in the sample as h, the number of mutation groups as m, and the representation of Y as haplotypes and mutation groups as Y h×m .

2.3 Representing Y h×m as a gene tree Y h×m (Figure 2A) can alternatively be represented as a gene tree or perfect phylogeny [START_REF] Gusfield | Efficient algorithms for inferring evolutionary trees[END_REF][START_REF] Griffiths | Sampling theory for neutral alleles in a varying environment[END_REF]. This representation relies on our assumption of the infinite sites mutation model in which, if a site mutates once in a given lineage, all descendants of that lineage also have the mutation and no other individuals carry that mutation. The gene tree is a graphical representation of the haplotypes (as tips) arranged by their patterns of shared mutations. The haplotype data summarized in Figure 2A corresponds to the gene tree given in Figure 2B. Details of the correspondence between haplotype data and gene tree are listed below, and an additional example is given in Figure 13 (Appendix E).

A gene tree for a matrix Y h×m of h haplotypes and m mutation groups is a rooted tree T with h leaves and at least m edges, such that (Figure 2B):

1. Each row of Y h×m corresponds to exactly one leaf of T . The black numbers at leaf nodes in Figure 2B are the haplotype frequencies.

2. Each mutation group of Y h×m is represented by exactly one edge of T , which is labeled accordingly (letters in Figures 2A and2B). The red numbers along edges in Figure 2B give the cardinality of each mutation group (i.e. the number of segregating sites in this group; see Figure 2A). Some external edges (edges subtending leaves) may not be labeled, indicating that they do not carry additional mutations to their parent edge. This happens when the other edges emanating from its parent node necessarily correspond to other mutation groups.

3. Edges are placed in the gene tree in such a way that each path from the root to a leaf fully describes a haplotype. Edges corresponding to shared mutations between several haplotypes are closest to the root. For example, in Figure 2B, haplotype 4 corresponds to the leaf at which one arrives starting from the root and going along edges a, g and f; in contrast, haplotype 7 corresponds to the leaf at which one arrives going from the root along edge b. Thus, the labels and the numbers associated with the edges along the unique path from the root to a leaf exactly specify a row of Y h×m . Dan Gusfield's perfect phylogeny algorithm [START_REF] Gusfield | Efficient algorithms for inferring evolutionary trees[END_REF] transforms the sequence data Y h×m into a gene tree and this transformation is one-to-one. We note that the perfect phylogeny T or gene tree is not the same as the genealogy g. While a genealogy is a bifurcating tree of individuals of the sample, the gene tree is a multifurcating tree of haplotypes.

Tajima's genealogies

Our method of computing the probability of the recoded data, Y h×m , uses ranked tree shapes rather than fully labeled histories. We refer to these ranked tree shapes as Tajima's genealogies but note they have also been called unlabeled rooted trees [START_REF] Griffiths | Unrooted genealogical tree probabilities in the infinitely-manysites model[END_REF] and evolutionary relationships [START_REF] Tajima | Evolutionary relationship of DNA sequences in finite populations[END_REF]. In Tajima's genealogies, only the internal nodes are labeled and they are labeled by their order in time. Tajima's genealogies encode the minimum information needed to compute the probability of data, Y h×m which consists of nested sets of mutations, without any approximations. In Figure 1A for example, it matters only that mutation group "e" occurs on a subgroup of the individuals who carry a mutation group "a" and that this is different than the subgroups carrying "c", "d" and "f". No other labels matter because individuals are exchangeable in the population model we assume.

D DAG representation Z 0 Z 1 Z 2 Z 3 Z 4 Z 5 Z 6 Z 7 Z 8 Z 9 Z 10 Z 11 B Perfect Phylogeny 2 2 2 2 2 2 2 1 1 2 2 2 3 1 1 1 1 1 1 1 Number of mutations Frequency Z j =      (D j , X j , A j ) j ∈ V I (D j , A j ) j = 0 (D j , X j ) j ∈ V L !"#$
This represents a dramatic coarsening of tree space compared to the classical leaf-labeled binary trees of Kingman's coalescent. The number of possible ranked tree shapes for a sample of size n corresponds to the n-th term of the sequence A000111 of Euler zig-zag numbers [START_REF] Disanto | Exact enumeration of cherries and pitchforks in ranked trees under the coalescent model[END_REF] whereas the number of labeled binary tree topologies is n!(n -1)!/2 n-1 . As can be seen from Figure 1B, this provides a much more efficient way to integrate over the key hidden variable, the unknown gene genealogy of the sample, when computing likelihoods.

We model this hidden variable using the vintaged and sized coalescent [START_REF] Sainudiin | Finding the best resolution for the kingman-tajima coalescent: theory and applications[END_REF] which corresponds exactly to this coarsening of Kingman's coalescent. As can be seen in Figure 1A, we assign vintages/labels 2 through n starting at the root of the tree and moving toward the present, so that the node created by the final splitting event, which is also the first coalescence event looking back in the ancestry of the sample, is labeled n. We write t k for the time of node k, measured from the present back into the past. We set t n+1 := 0 to be the present time. Then during the interval [t k+1 , t k ) the sample has exactly k extant ancestors, for k ∈ {2, . . . , n}.

The coarsening of the tree topology does not change the law of the times between two coalescence events. Thus, conditional on the effective population size trajectory (N (t)) t≥0 and the time t k+1 at which the number of ancestors to the sample decreases to k, the distribution of the time during which the sample has k ancestors is given by [START_REF] Slatkin | Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations[END_REF], where

Pr t k -t k+1 |t k+1 , (N (t)) t≥0 = C k N (t k ) exp - t k t k+1 C k N (t) dt (2)
C k = k 2 .
Writing the density at t = (t 2 , t 3 , . . . , t n ) of the vector of coalescence times as a product of conditional densities, we obtain

Pr t | (N (t)) t≥0 = n k=2 Pr t k -t k+1 | t k+1 , (N (t)) t≥0 . (3)
We use a lower-triangular matrix denoted F to represent Tajima's genealogies; see Appendix A. The probability of a ranked tree shape was derived independently in [START_REF] Sainudiin | Finding the best resolution for the kingman-tajima coalescent: theory and applications[END_REF] and [START_REF] Palacios | Bayesian nonparametric inference of population size changes from sequential genealogies[END_REF]. Specifically, for every ranked tree shape F with n leaves,

Pr(F) = 2 n-c-1 (n -1)! , ( 4
)
where c is the number of cherries in F (i.e., nodes subtending two leaves; c = 3 in Figure 10A).

Note that this probability is independent of the effective population size trajectory since the choice of the pair of lineages that coalesce during an event is independent of (N (t)) t≥0 (recall that in Kingman's coalescent, the coalescing pair is chosen uniformly at random among all possible pairs).

Since the distribution of Tajima's genealogies g T = (F, t) conditional on (N (t)) t≥0 can be factored as the product of the probability of the ranked tree shape F and the coalescent times density, we arrive at

Pr(g T | (N (t)) t≥0 ) = 2 n-c-1 (n -1)! n k=2 C k N (t k ) exp - t k t k+1 C k N (t) dt .
(5)

An augmented data representation using directed acyclic graphs

A key component of BESTT is the calculation of the conditional likelihood Pr(Y|g T , µ). We compute the conditional likelihood recursively over a directed acyclic graph (DAG) D. Our DAG exploits the gene tree representation T of the data (Figure 2B), incorporates the branch length information of the Tajima's genealogy g T (Figure 2C) and facilitates the recursive allocation of mutations to the branches of g T . Here we detail the construction of the DAG.

We construct the DAG using three pieces of information: the observed gene tree T , a given Tajima's genealogy g T and a latent "allocation" of mutations along the branches of the Tajima's genealogy (Figure 3). An allocation refers to a possible mapping (compatible with the data)

of the observed numbers of mutations (red numbers in Figure 2B) to branches in the Tajima's genealogy. Figure 3A shows one possible mapping for the Tajima's genealogy in Figure 2C; usually this mapping is not unique. Our construction of D enables an efficient recursive consideration of all possible allocations of mutations along g T when computing the conditional likelihood Pr(Y | g T , µ).

Constructing the DAG D. The graph structure of our DAG D = {Z, E} (Figure 2D) with nodes Z and edges E is constructed from a gene tree T . The number of internal nodes in the DAG D is the same as the number of internal nodes in T . However, sister leaf nodes in T with the same number of descendants are grouped together in D and leaf nodes descending from edges with no mutations are treated as singletons grouped together in D. For example, the leaves in Figure 2B subtending from edges i and j are grouped into Z 6 in Figure 2D, as they both have haplotype frequency 2. However, the leaves subtending from the e and f edges are not grouped (and correspond to Z 8 and Z 9 in the DAG Figure 2D) since they have respective haplotype frequencies 2 and 1.

We label the root node of D as Z 0 and increase the index i of each node Z i from top to bottom, moving left to right. For i < j, we assign a directed edge E i,j if the node in T corresponding to Z i is connected to the node in T corresponding to Z j . The index set of internal nodes in D is denoted by V I and the index set of leaf nodes is denoted by V L .

Information carried by the nodes in D. Each node in D represents a vector, Z j , which includes number of descendants, number of mutations and latent allocation of mutations. Although the number of descendants and number of mutations are part of the observed data, the allocation of mutations can be seen as a random variable, for ease of exposition, we use capital letters to denote all three types of information. We define the vector Z j as follows:

Z j =      (D j , X j , A j ) j ∈ V I , (D j , A j ) j = 0 (the root node), (D j , X j ) j ∈ V L ,
where D j denotes the number of descendants of (i.e., of sampled sequences subtended by) node Z j , X j denotes the number of mutations separating Z j from its parent node, and A j denotes the allocation of mutations along g T (described in detail below). The number of descendants D j is thus the number of individuals/sequences descending from node Z j (this information is part of T ). For internal nodes, X j records the cardinality of a mutation group, represented as a red number along the edge E i,j of T in Figure 2B, where i is the index of the parent node of Z j . Leaf nodes in D may correspond to more than one leaf nodes in T , namely any sister nodes with the same number of descendants. In this case, X j is a vector with the cardinalities of the corresponding mutation groups (see for example node Z 6 in Figure 3B). In order to keep the DAG construction simple, we only allow groupings of leaf nodes and not of internal nodes with identical descendants carrying identical numbers of mutations. We note that, in principle, it would be possible to compress the number of internal nodes of the DAG by exploiting all the symmetries observed in the data.

Allocation of mutation groups along g T . The latent allocation variables {A j } determine a possible correspondence between subtrees in g T and nodes in D: in particular, A j indicates the branches in g T that subtend the subtrees corresponding to nodes When leaf nodes represent more than one trees, such as Z 4 in Figure 4B, the random variable X j is the vector X j = (X j,1 , X j,2 , . . . , X j,sj ) that denotes the number of mutations along the branches that subtends from the tree node j that have D j descendants, and s j is the number of edges subtending from Z j .

{Z k } if {Z k } are child nodes of Z j . B DAG corresponding to A Z j =      (D j , X j , A j ) j ∈ V I (D j , A j ) j = 0 (D j , X j ) j ∈ V L 16
z 3 = (d 3 = 2, x 3 = 2) z 4 = (d 4 = 2, x 4 = (1, 1)) z 5 = (d 5 = 3, x 5 = 1, a 5 = (b 16 , b 24 )) z 6 = (d 6 = 2, x 6 = (1, 2)) z 7 = (d 7 = 3, x 7 = 1, a 7 = (b 15 , b 25 )) z 8 = (d 8 = 2, x 8 = 2) z 9 = (d 9 = 1, x 9 = 3) z 10 = (d 10 = 2, x 10 = 1) z 11 = (d 11 = 1, x 11 = 1) Z 0 Z 1 Z 2 Z 3 Z 4 Z 5 Z 6 Z 7 Z 8 Z 9 Z 10 Z 11 b 24 b 25 z 1 = (d 1 = 7, x 1 = 1, a 1 = ((b 12 , b 9 ), b 10 )) z 2 = (d 2 = 7, x 2 = 1, a 2 = ((b 8 , b 13 ), b 11 ))
Allocations of mutations to branches are usually not unique and computation of the conditional likelihood Pr(Y | g T , µ) requires summing over all possible allocations. In Figure 3A we show one such possible allocation of the mutation groups of the gene tree in Figure 2B along the Tajima's genealogy in Figure 2C. For example, mutation group "a" in Figure 2B with cardinality 1 (number in red) is a mutation observed in 7 individuals (sum of black numbers of leaves descending from edge marked a). This same mutation group, "a", is shown as a red number 1 in Figure 3A allocated to branch b 5 . If Z j is an internal node, the number of mutations X j is denoted as a vector of length 1.

If Z j is a leaf node, X j can be a vector of length greater than 1. Details on notation for allocations can be found in Appendix B.

Computing the conditional likelihood

Under the infinite-sites mutation model, mutations are superimposed independently on the branches of g T as a Poisson process with rate µ. In order to compute Pr(Y | g T , µ) = Pr(T | g T , µ) we marginalize over the latent allocation information in the directed acyclic graph D; that is, we sum over all possible mappings of mutations in T to branches in g T as follows:

Pr(Y | g T , µ) = A 0 A 1 .
. .

An I Pr(D | g T , µ) = A 0 A 1 .
. .

An I Pr(Z 0 , . . . , Z n I +n L | g T , µ) = A 0 A 1 .
. .

An I n I +n L i=1 Pr(Z i | Z pa(i) , g T , µ)
where

n I = |V I |, n L = |V L |, pa(i)
denotes the index of the parent of node i in D and we set P (Z 0 | g T , µ) = 1 because it is assumed that there are no mutations above the root node and the length of the root branch l 2 = 0. Writing L for the tree length of g T (i.e., the sum of the lengths of all branches of g T ) and factoring out a global factor e -µL (due to the Poisson distribution of mutations across the genealogy) from each of the above products over i ∈ {1, . . . , n I + n L }, we have

Pr(Z i = z i | z pa(i) , g T , µ) =      Pr(X i = x i | a pa(i) = b j , g T , µ) ∝ (µl j ) x i if |x i | = 1, Pr X i = (x i1 , . . . , x ik ) | a pa(i) = (b j 1 , . . . , b j k ), g T , µ ∝ s∈Π(x i ,k) k m=1 (µl jm ) sm if |x i | = k > 1,
where Π(x i , k) is the set of all permutations of x i = {x i1 , . . . , x ik } divided into m i groups of different sizes. The number of different permutations of the k values of x i divided into m i groups of sizes

k 1 , . . . , k m i is |Π(x i , k)| = k! m i j=1 k j ! (6)
For example, assume that x i = {2, 2, 2, 0, 3, 3} and a pa

(i) = (b 3 , b 4 , b 5 , b 6 , b 7 , b 8
) with branch lengths {l 3 , l 4 , l 5 , l 6 , l 7 , l 8 }. In this case, k 1 = 3 because there will be 3 branches with 2 mutations, k 2 = 1 because there will be 1 branch with 0 mutations and k 3 = 2 because there will be 2 branches with 3 mutations. The number of permutations of k = 6 mutations groups divided into m i = 3 groups with cardinalities 2, 0,

3 of sizes 3, 1, 2 is 6!/(3!1!2!) = 60.
The conditional likelihood Pr(Y | g T , µ) is calculated via a backtracking algorithm (Appendix C). The algorithm marginalizes the allocations by traversing the DAG from the tips to the root. The pseudocode and an example can be found in the Appendix C.

The case of unknown ancestral states

Up to now, we have assumed that the ancestral state was known at every segregating site. The representation of the data Y that we use in this case records the cardinalities of each mutation group and the genealogical relations between these groups, but does not assign labels to the sequences.

Hence, in the terminology of [START_REF] Griffiths | Unrooted genealogical tree probabilities in the infinitely-manysites model[END_REF], our data corresponds to an unlabeled rooted gene tree.

When the ancestral types are not known, the data (now denoted Y 0 ) may be represented as an unlabeled unrooted gene tree. By the remark following Equation ( 1) in [START_REF] Griffiths | Unrooted genealogical tree probabilities in the infinitely-manysites model[END_REF], if s is the number of segregating sites, then there are at most s + 1 unlabeled rooted gene trees that correspond to the unrooted gene tree of the observed data (R(Y 0 )). By the law of total probability (see also Equation ( 10) in [START_REF] Griffiths | Unrooted genealogical tree probabilities in the infinitely-manysites model[END_REF]), the conditional likelihood of Y 0 can be written as the sum over all compatible unlabeled rooted gene trees Y (i) of the probability of Y (i) conditionally on g T . That is:

Pr Y 0 | g T , µ = R(Y 0 ) i=1 P Y (i) | g T , µ , (7)
where each of the Y (i) corresponds to a unique unlabeled rooted gene tree compatible with the unrooted gene tree Y 0 and R(Y 0 ) denotes the number of those unlabeled rooted gene trees. In the following sections, we shall assume that the ancestral type at each site is known.

Bayesian inference of the effective population size trajectory

Our posterior distribution of interest is

Pr(γ, g T , τ | Y, µ) ∝ Pr(Y | g T , µ)Pr(g T | γ)Pr(γ | τ )Pr(τ ), (8)
where (log N (t)) t≥0 = (γ(t)) t≥0 ∼ GP(0, C(τ )) has a Gaussian process prior with mean 0 and covariance function C(τ ) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. This specification ensures (N (t)) t>0 is non-negative. In our implementation, we assume a regular geometric random walk prior, that is,

γ 1 = log N (t * 1 ), . . . , γ B = log N (t * B ) at B regularly spaced time points in [0, T ] with Cov[γ i , γ j ] = Cov log N (t * i ), log N (t * j ) = τ min(t * i , t * j ).
The parameter τ is a length scale parameter that controls the degree of regularity of the random walk. We place a Gamma prior with parameters α = .01 and β = .001 on τ , reflecting our lack of prior information in terms of high variance about the smoothness of the logarithm of the effective population size trajectory.

We approximate the posterior distribution of model parameters via a MCMC sampling scheme.

Model parameters are sampled in blocks within a random scan Metropolis-within-Gibbs framework.

Our algorithm initializes with the corresponding Tajima genealogy of the UPGMA estimated tree implemented in phangorn [START_REF] Schliep | phangorn: phylogenetic analysis in R[END_REF]. Given an initial genealogy, our algorithm initializes N e and τ with the method of (Palacios and Minin, 2012) implemented in phylodyn [START_REF] Karcher | phylodyn: an r package for phylodynamic simulation and inference[END_REF]. We then proceed to generate (1) a sample of the vector of effective population sizes and precision parameter as described in section 2.8.2, (2) a sample of the vector of coalescent times as described in 2.8.3 and 2.8.4 where we modify a single coalescent time, and (3) a sample of ranked tree shape as described in 2.8.1 in each iteration. To summarize the effective population size trajectory, we compute the posterior median and 95% credible intervals pointwise at each grid point in [0, T ], were T is the maximum time to the most recent common ancestor sampled.

1 2 3 k k-1 1 2 k 3 4 k-1 2 3 1 k k-1 1 3 2 k k-1 1 2 k 3 4 k-1

Metropolis-Hastings updates for ranked tree shapes

There is a large literature on local transition proposal distributions for Kingman's topologies [START_REF] Kuhner | Maximum likelihood estimation of population growth rates based on the coalescent[END_REF][START_REF] Rannala | Bayes estimation of species divergence times and ancestral population sizes using dna sequences from multiple loci[END_REF][START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF][START_REF] Whidden | Quantifying mcmc exploration of phylogenetic tree space[END_REF][START_REF] Aberer | An efficient independence sampler for updating branches in bayesian markov chain monte carlo sampling of phylogenetic trees[END_REF]. In this paper, we adapted the local transition proposal of [START_REF] Markovtsova | The age of a unique event polymorphism[END_REF] to Tajima's topologies. We briefly describe the scheme below and provide a pseudocode algorithm in Appendix C (Algorithm 3).

Given the current state of the chain {γ, τ, g T } = {γ, τ, F n , t}, we propose a new ranked tree shape F * in two steps: [START_REF]Z 8 -Z 10 -Z 3 -Z 6 -Z 4 -Z 7 -Z 5 -Z 4 -Z 6 -Z 1 -Z 2 -Z 1 -Z 2 -Z 0 2. Z 8 -Z[END_REF] we first sample a coalescent interval e k = (t k+1 , t k ) uniformly at random, where k ∼ U({3, . . . , n}). Note that we will never select the interval (t 3 , t 2 ) at the top of the tree (see Figure 10A). Given k, we focus solely on the coalescent events at times t k and t k-1 . For step (2), there are two possible scenarios. Case A: The lineage created at time t k , labeled k, coalesces at time t k-1 (first row of Figure 4A). Case B: Lineage k does not coalesce at time t k-1 (Figure 4B).

In Case A, we choose a new pair of lineages at random to coalesce at time t k from the 3 lineages subtending k and k -1 (excluding k), and we coalesce the remaining lineage with k at t k-1 (F * n,1 and F * n,2 in Figure 4). In Case B, we invert the order of the coalescent events; that is, the two lineages descending from k are set to coalesce at time t k-1 and lineages descending from k -1 are set to coalesce at time t k . (F * n,3 a in Figure 4). Note that the numerical labels 1, 2, 3 are included to clarify the picture: lineages subtending both Case A and Case B can be either labeled (if there is a vintage subtending that lineage) or not (if there is a singleton). The transition probability

q(F * n | F n ) is
given by the product of the probabilities of the two steps. The new ranked tree shape F * n is accepted with probability given by the Metropolis-Hastings ratio defined below:

a Fn = min 1 , Pr(Y | F * n , t, µ)Pr(F * n )q(F n | F * n ) Pr(Y | F n , t, µ)Pr(F n )q(F * n | F n ) (9)
We note that our proposal can result in the same ranked tree shape. However, we tested alternative proposals that precluded this event and we did not find any notable difference in the overall performance of the MCMC algorithm.

Split Hamiltonian Monte Carlo updates of (γ, τ )

To make efficient joint proposals of γ and τ , we use the Split Hamiltonian Monte Carlo method proposed by [START_REF] Lan | An efficient Bayesian inference framework for coalescent-based nonparametric phylodynamics[END_REF]. Conditioned on g T , the target density becomes π(γ, τ ) ∝ Pr(t | γ)Pr(γ | τ )Pr(τ ). This is the same target density implemented in [START_REF] Karcher | phylodyn: an r package for phylodynamic simulation and inference[END_REF] for fixed coalescent times t.

Hamiltonian Monte Carlo updates of coalescent times

Given the current state {γ, τ, g T } = {γ, F n , t, τ }, we propose a new vector of coalescent times with target density π(t ) ∝ P (Y | F n , t , µ)P (t | γ) by numerically simulating a Hamilton system with Hamiltonian

H(log(t ), s) = -log π(log(t )) + 1 2 s T Ms, ( 10
)
where s is the momentum vector assumed to be normally distributed. The system evolves according to:

∂s ∂x = ∇ log π(log(t ))) ∂t ∂x = Ms (11)
We use the leapfrog method [START_REF] Neal | MCMC using Hamiltonian dynamics[END_REF] with step size and a p Poisson with mean 10 distributed number of steps to simulate the dynamics from time x = 0 to x = p . Each leapfrog step of size follows the trajectory:

s x+ /2 = s x + 2 ∇ log π(log(t x ))) t x+ = t x + Ms x+ /2 s x+ = s x+ /2 + 2 ∇ log π(log(t x+ ))) (12)
For our implementation, we set the mass matrix M = I, the identity matrix. We simulate the Hamiltonian dynamics of the logarithm of times to avoid proposals with negative values. Solving the equations of the Hamilton system requires calculating the gradient of the logarithm of the target density with respect to the vector of log coalescent times. The gradient of the log conditional likelihood (score function) is calculated at every marginalization step in the algorithm for the likelihood calculation.

At the beginning of Section 2.8, we described how we assume a regular geometric random walk prior on (N (t)) t≥0 at B regularly spaced time points in [0, T ]. Ideally, the window size T must be at least t 2 , the time to the most recent common ancestor (TMRCA). However, t 2 is not known.

Our initial values of coalescent times t are obtained from the UPGMA implementation in phangorn [START_REF] Schliep | phangorn: phylogenetic analysis in R[END_REF] with times properly rescaled by the mutation rate, and we set T = t 2 . We initially discretize the time interval [0, T ] into B intervals of length T /(B -1). As we generate new samples of t, we expand or contract our grid according to the current value of t 2 by keeping the grid interval length fixed to T /(B -1), effectively increasing or decreasing the dimension of γ.

Local updates of coalescent times

In addition to HMC updates of coalescent times, we propose a move of a single coalescent time (excluding the TMRCA t 2 ) chosen uniformly at random and sampled uniformly in the intercoalescent 4. A constant trajectory:

N (t) = 1
Given a genealogy of length L = n j=2 j(t jt j+1 ), where t jt j+1 is the intercoalescent length while there are j lineages, we drew the total number of mutations (segregating sites) s according to a Poisson distribution with parameter µL. We then placed the mutations uniformly at random along the branches of the genealogy. For each of the s mutations, we assigned the mutant type to individuals descending from the branch where the mutation occurred and the ancestral type otherwise.

We summarize our posterior inference N (t) by the posterior median and 95% Bayesian credible intervals after 200 thousand iterations and thinned every 10 iterations with 100 iterations of burn in. Our initial number of change points for N (t) was set to 50 over the time interval (0, t 2 ), where t 2 is the initialized time to the most recent common ancestor; however, over the course of MCMC iterations, this number could increase or decrease according to the posterior distribution of t 2 .

We assess accuracy and precision of our estimates using the sum of relative errors (SRE)

SRE = k i=1 | N (ω i ) -N (ω i )| N (ω i ) , ( 17
)
where N (ω i ) is the estimated effective population size trajectory at time ω i . Second, we computed the mean relative width as

M RW = k i=1 | Nup (ω i ) -Nlo (ω i )| kN (ω i ) , (18)
where Nup (ω i ) corresponds to the 97.5% upper limit and Nlo (ω i ) corresponds to the 2.5% lower limit of the estimated posterior distribution of N (ω i ). In addition, we measured how well the 95% credible intervals cover the truth and compute the envelope measure, EN V :

EN V = k i=1 1( Nlo (ω i ) ≤ N (ω i ) ≤ Nup (ω i )) k (19)
We first simulated 3 datasets of n = 10 individuals with an average number of 100 segregating sites under different types of population size trajectories: constant, exponential growth and instantaneous growth. Results are depicted in the first column of Figure 8. Posterior medians and 95% credible intervals are shown as black curves and gray shaded areas respectively. The trajectory used to simulate the data is depicted as a dashed line. Figure 8 shows that our BESTT method recovers the constant and exponential growth trajectories very well but the instantaneous growth scenario is less accurate and with high uncertainty (wide credible intervals). In all three cases, our envelope measure is above 95%. Performance measures on all simulations are summarized in Table 1.

We analyzed the effect of increasing the number of segregating sites, the number of samples and the number of independent genealogies on posterior inference with BESTT. In all three cases, we expect our method to better recover the truth. Figure 5 shows our results on simulated data under a population size trajectory with instantaneous growth (Equation 16) of n = 10 individuals with 31, 63 and 120 segregating sites. As expected, our method recovers the truth with higher precision (MRW) and accuracy (SRE) when we increase the number of segregating sites. Increasing the number of segregating sites may result in more constraints in the gene tree. For n = 10, there are 7936 possible ranked tree shapes, however for the datasets simulated with 31, 63 and 102 segregating sites, there are only 2582 ± 32, 2670 ± 34 and 556 ± 7 ranked tree shapes compatible with their sequences under a population size trajectory with instantaneous growth (dashed lines). s is the number of segregating sites. Posterior medians are depicted as solid black lines and 95% Bayesian credible intervals are depicted by shaded areas. corresponding gene trees. These numbers were estimated by importance sampling [START_REF] Cappello | Sequential importance sampling for multi-resolution kingman-tajima coalescent counting[END_REF]. As another performance assessment, we simulated datasets from a population size trajectory with instantaneous growth with varying number of samples. We simulated datasets with n = 10, 25 and 35 samples with 215 expected number of segregating sites. Our results depicted in Figure 6 show that our method performs better in terms of SRE and MRE when the number of samples increases. Similarly, precision (MRW) and accuracy (SRE) increases when inference is done from a larger number of independent datasets. Finally, Figure 7 shows our results from 1, 5 and 10 datasets simulated from 1, 5 and 10 independent genealogies of 10 individuals with a population size trajectory of growth followed by a constant period (Equation 15). As expected, our method's performance substantially increases by increasing the number of independent datasets.

Comparison to other methods

To our knowledge, there is no other method for inferring (variable) effective population size over time from haplotype data that assumes the infinite sites mutation and a nonparametric prior on N (t), therefore we cannot have a direct comparison of our method to others. Moreover, our method is the only one that explicitly averages over Tajima genealogies instead of Kingman genealogies.

BEAST [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF]) is a program for analyzing molecular sequences that uses MCMC to average over the Kingman tree space and it is therefore a good reference for comparison to our method. We compared our results to the Extended Bayesian Skyline Plot method (EBSP) [START_REF] Heled | Bayesian inference of population size history from multiple loci[END_REF] and the Skygrid method [START_REF] Gill | Improving bayesian population dynamics inference: A coalescent-based model for multiple loci[END_REF] implemented in BEAST.

Since the infinite sites mutation model is not implemented in BEAST, we first converted our simulated sequences of 0s and 1s to sequences of nucleotides by sampling s ancestral nucleotides uniformly on {A, T, C, G} and assigning one of the remaining 3 types uniformly at random to be the mutant type. This corresponds to a simulation of the Jukes-Cantor mutation model [START_REF] Jukes | Evolution of protein molecules[END_REF] that is currently implemented in BEAST.

We compare the results of BESTT to those of BEAST EBSP and Skygrid [START_REF] Drummond | Bayesian coalescent inference of past population dynamics from molecular sequences[END_REF][START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF] in Figure 8. We note that results from BEAST are generated from 10 million iterations and thinned every 1000 iterations, while results from BESTT are generated from 200 thousand iterations.

We compared our point estimates N (t) from all methods to the ground truth for each simulation (Table 2). In two cases, BESTT has better envelope than BEAST. For the exponential growth simulation (Figure 8, second row) the BEAST EBSP result has better SRE and MRW, however, the credible intervals are uneven with very wide intervals at the ends. In all cases, the BEAST Skygrid results have wider credible intervals. For the instantaneous growth simulation (Figure 8, third row), BEAST EBSP did not generate many simulations beyond the time point 0.06, for this reason we recomputed the performance statistics for the overlapping time interval (0, 0.06). In this interval, BESTT outperforms both methods implemented in BEAST in terms of envelope and SRE. The last column of Figure 8 shows the posterior distribution of the time to the most recent common ancestor (TMRCA). For the case of constant population size, the true value of the TMRCA is contained in the 95% BCI estimated with BESTT but it is not contained in the 95% BCIs estimated with the two methods implemented in BEAST. In the exponential growth simulation, the true TMRCA is contained in the 95% BCIs estimated with the three methods and the instant growth method, the true TMRCA is not contained in the 95% BCIs of the three methods. We note that BEAST Bayesian Skygrid [START_REF] Gill | Improving bayesian population dynamics inference: A coalescent-based model for multiple loci[END_REF] is a more comparable method to BESTT since it assumes Gaussian process priors on log N (t) like BESTT.

Computational performance of BESTT

BESTT approximates the posterior distribution (a) Pr((N (t)) t≥0 , g T , τ | Y, µ), where g T is a

Tajima's genealogy instead of (b) Pr((N (t)) t≥0 , g, τ | Y * , µ), where g is a Kingman's genealogy [START_REF] Cappello | Sequential importance sampling for multi-resolution kingman-tajima coalescent counting[END_REF]. That is, we are required to estimate the posterior of a smaller number of trees. For this reason, we argue that Tajima's coalescent is a more efficient model than Kingman's coalescent for estimating the posterior distribution of (N (t)) t≥0 . However, a single conditional likelihood calculation Pr(Y | g T , µ) requires the sum over all possible allocation of mutation groups to branches of g T . Our algorithm only accounts for allocations constrained by the DAG and the ranked tree shape of g T . For the data depicted in Figure 2A,B and g T of Figure 2C, there are only 8 different possible allocation "paths" of all mutation groups to branches. In Appendix C we detail how our algorithm finds these paths.

The number of paths depends on the number of subtrees with the same family size path in the DAG and in the ranked tree shape. In the best case, our algorithm will find a path in O(no),

where no is the number of nodes in the gene tree. In general, the number of allocation paths will be much smaller than the number of labeled trees compatible with a ranked tree shape. In our implementation, we estimate posterior (a) with MCMC. The main difference between our MCMC algorithm and the MCMC algorithm implemented in BEAST is the tree topology sampler. While our MCMC algorithm explores the space of ranked tree shapes with local move proposals of ranked tree shapes, BEAST explores the space of labeled, ranked tree shapes with local move proposals of labeled trees. A formal assessment of the efficiency of our MCMC algorithm and its comparison to the MCMC implementation in BEAST is beyond the scope of this manuscript and subject of future research.

Inferring human population demography from mtDNA

We selected n = 35 samples of mtDNA at random from 107 Yoruban individuals available from the 1000 Genomes Project phase 3 (The 1000 Genomes Project Consortium, 2015). We retained the coding region: base pairs 576 -16, 024 according to the rCRS reference of Human Mitochondrial DNA [START_REF] Anderson | Sequence and organization of the human mitochondrial genome[END_REF][START_REF] Andrews | Reanalysis and revision of the Cambridge reference sequence for human mitochondrial dna[END_REF] and removed 38 indels. Of the 260 polymorphic sites, we retained 240 sites compatible with the infinite sites mutation model. The final file is available in https://github.com/JuliaPalacios/phylodyn. To encode our data as 0s and 1s, we use the inferred root sequence RSRS of [START_REF] Behar | A "copernican" reassessment of the human mitochondrial dna tree from its root[END_REF] to define the ancestral type at each site. To rescale our results in units of years, we assumed a mutation rate per site per year of 1.3 × 10 -8 (Rebolledo-Jaramillo et al., 2014). We compare our results with the Extended Bayesian Skyline method [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF] implemented in BEAST in Figure 9. When applying BEAST, we assumed the Jukes-Cantor mutation model. Both methods detect an inflection point around 20kya followed by exponential growth. The mean time to the most recent ancestor (TM-RCA) inferred for these YRI mtDNA samples with BESTT is around 170kya with a 95% BCI of (142868, 207455), while the mean TMRCA inferred with BEAST is around 160kya with a 95% BCI of (133239,196900). In Appendix D, we include two more comparisons of BESTT and BEAST.

Discussion

The size of emergent sequencing datasets prohibits the use of standard coalescent modeling for inferring evolutionary parameters. The main computational bottleneck of coalescent-based inference of evolutionary histories lies in the large cardinality of the hidden state space of genealogies. In the standard Kingman coalescent, a genealogy is a random labeled bifurcating tree that models the set of ancestral relationships of the samples. The genealogy accounts for the correlated structure induced by the shared past history of organisms and explicit modeling of genealogies is fundamental for learning about the past history of organisms. However, the genomic era is producing large datasets that require more efficient approaches that efficiently integrate over the hidden state space of genealogies.

In this manuscript we show that a lower resolution coalescent model on genealogies, the "Tajima's coalescent", can be used as an alternative to the standard Kingman coalescent model. In particular, we show that the Tajima coalescent model provides a feasible alternative that integrates over a smaller state space than the standard Kingman model. The main advantage in Tajima's coalescent is to model the ranked tree topology as opposed to the fully labeled tree topology as in Kingman's coalescent.

A priori, the cardinality of the state space of ranked tree shapes is much smaller than the cardinality of the state space of labeled trees. However, in this manuscript we show that when the Tajima coalescent model is coupled with the infinite sites mutation model, the space of ranked tree shapes is constrained by the data and the reduction on the cardinality of the hidden state space of Tajima's trees is even more pronounced than expected.

supported by the National Science Foundation CAREER Award DBI-1452622 to S.R. 0 1 3 0 1 2 4 0 1 2 3 5 0 1 2 3 4 6 0 1 2 3 3 5 7 0 1 1 2 2 4 6 1A) with coalescent events ranked from 2 at time t 2 to n at time t n . Right: The corresponding F n matrix, with n = 8, that encodes the ranked tree shape information of the Tajima's genealogy on the left. F i,j denotes the number of lineages that do not coalesce in the time interval (t i+1 , t j ). In particular, F n,i for i in {2, 3, . . . , n} denotes the number of singletons (external branches that have not coalesced) in the time interval (t i+1 , t i ).

F =             0 0 2
Algorithm 1 Calculate Likelihood(F path, times, GeneT ree) procedure Input: GeneT ree , F P ath Output: Log Likelihood LL 1: Initiate pool to be the set of leave nodes of GeneT ree with at least one descendant. Initiate LL and index to be zero. Initiate current path to be empty. 2: Call CalcLL recursive(LL, index, current path, F path, times, Genetree). for node in tree do 3:

Calculate log likelihood based on times and number of mutations of node in current path.

4:

Accumulate to total log likelihood LL Check compatibility of the node, according to the given F path.

9:

if node is compatible with F path then 10:

Update node by assigning it to the current step in F path 11:

Update pool. If a node has been mapped entirely, remove node from pool, update its parent node, and potentially add parent node to pool if parent node has not been entirely assigned.

12:

Append this node to current path end for 17: end if To illustrate our algorithms 1 and 2, we use our example of Figures 2 and3. Algorithm 1 initiates the pool with nodes Z 3 , Z 4 , Z 6 , Z 8 , Z 10 . Then, Algorithm 2 cycles through this list. Assume the first node is Z 8 . This node has d 8 = 2 descendants and the ancestry is Z 8 -Z 5 -Z 1 -Z 0 with sizes: 2 -3 -7 -16. On the other hand, the first coalescent event (from present to past) labeled 16 in g T (Figure 3A) has ancestry with sizes: 2 -3 -7 -16. Therefore, this node is compatible. The node is removed from the pool, its parent node added to the pool and Z 8 is assigned to the path. At this time current path = Z 8 and the pool becomes: Z 3 , Z 4 , Z 5 , Z 6 , Z 10 . The algorithm then cycles through this list and picks Z 10 . This node has size ancestry 2 -3 -7 -16. On the other hand the second coalescent event labeled 15 has size ancestry: 2 -3 -5 -7 -9 -16. Since the node size ancestry is contained in the second coalescent event's size ancestry, this node is compatible. The current path becomes current path = Z 8 -Z 10 and pool becomes Z 3 , Z 4 , Z 5 , Z 6 , Z 7 . We continue this procedure until we reach the path current path = Z 8 -Z 10 -Z 3 -Z 6 -Z 4 -Z 7 -Z 5 -Z 4 -

Z 6 -Z 1 -Z 2 -Z 1 -Z 2 -Z 0 .
Once a path is found, the algorithm back tracks the path until there is one compatible node and the path continues to grow. A sequence of back tracking and growing is the following: to the root shows that this haplotype has a unique mutation corresponding to mutation group a.

We note that mutation group labels carry no information. We incorporate the labels in the Figure for ease of exposition. Since mutation group h has now multiplicity 2, the branch labeled h has now a red label 2.

Figure 2 :

 2 Figure 2: Data structures employed by our method, BESTT, for calculating the conditional likelihood of the data. A. Compressed data representation Y h×m of n = 16 sequences and s = 18 (columns, only the first 10 of which are shown), comprised of 9 haplotypes and 13 mutation groups. Rows correspond to haplotypes and each polymorphic site is labeled by its mutation group {a, b, c, ..., m}. B. Gene tree representation of the data in panel A. Red numbers indicate the cardinality of each mutation group (number of columns with the same label in panel A). Black letters indicate the mutation group (column labels in panel A), and black numbers indicate the frequency of the corresponding haplotype. C. A Tajima's genealogy compatible with the gene tree in panel B. Internal nodes are labeled according to order of coalescent events from the root to the tips. Coalescent event i happens at time t i and branches are labeled b i (see section 2.5 for details). D. A Directed Acyclic Graph (DAG) representation of the gene tree in panel B together with allocation of mutation groups along the branches of the Tajima's genealogy in panel C. V I denotes the set of internal nodes and V L the set of leaf nodes. A detailed description of the DAG is given in section 2.5.

  z 0 = (d 0 = 16, a 0 = (b 5 , b 4 , b 14 ))

b 5 1 1Figure 3 :

 13 Figure 3: DAG Construction. A. A Tajima's genealogy from Figure 2C with added allocation of mutations shown in red. B. The corresponding augmented DAG with allocation of mutations. At the root Z 0 , there are no mutations by convention. Node Z 0 has 16 descendants across 3 subtrees of 7, 7 and 2 descendants, corresponding to nodes Z 1 , Z 2 , Z 3 . These three subtrees subtend from b 5 , b 4 and b 14 , respectively, in g T (Figure 3A). Node Z 1 corresponds to the tree subtending from b 5 of size 7 with X 1 = 1 mutation along b 5 and subtends three subtrees from (b 12 , b 9 ) and b 10 . Subtrees subtending from (b 12 , b 9 ) are grouped together in leaf node Z 4 because they both have 2 descendants and have the same parent node.When leaf nodes represent more than one trees, such as Z 4 in Figure4B, the random variable X j is the vector X j = (X j,1 , X j,2 , . . . , X j,sj ) that denotes the number of mutations along the branches that subtends from the tree node j that have D j descendants, and s j is the number of edges subtending from Z j .

Figure 4 :

 4 Figure 4: Markov moves for topologies First row: possible coalescent patterns (Case A or Case B) for a given topology F n . Second row: possible Markov moves in Case A (F * n,1 and F * n,2 ) and Case B (F * n,3). k indexes the coalescent interval sampled. Numerical labels at the tips are added for convenience: conditionally on a given F n , tips can be labeled(vintage) or not (singleton). Figure is adapted from Figures 2, 3 and 4 of[START_REF] Markovtsova | The age of a unique event polymorphism[END_REF].

Figure 5 :

 5 Figure 5: Varying the number of segregating sites. Posterior inference from simulated data of n = 10

Figure 6 :

 6 Figure 6: Varying the number of samples under a population size trajectory with instantaneous growth. Posterior inference from simulated data of n = 10, 25 and 35 sequences under the population size trajectory with instantaneous growth. Shaded areas correspond to 95% credible intervals, solid lines to posterior median and dashed line to the truth.

Figure 7 :

 7 Figure7: Multiple independent datasets. Posterior inference from simulated data of n = 10 sequences under exponential followed by constant trajectory (eq. 15). (A) Inference from a single simulated dataset, (B) from 5 independently simulated datasets, and (C) from 10 independently simulated datasets. Shaded areas correspond to 95% credible intervals, solid black lines show posterior medians and dashed lines show the simulated truth.

Figure 9 :

 9 Figure 9: Posterior inference of female effective population size from 35 mtDNA samples from Yoruban individuals in the 1000 Genomes Project using our method BESTT (first plot) and the BEAST Extended Bayesian Skyline Plot (second plot). Posterior median curves are depicted as solid black lines and 95% credible intervals by shaded regions.

Figure 10 :

 10 Figure 10: Ranked tree shape Left: Example of a Tajima's genealogy (redrawn from Figure 1A) with

  3: return LL Algorithm 2 CalcLL recursive(LL, index, current path, F path, times, Genetree) procedure 1: if index = len(path) {When a complete path node is found} then 2:

13 :

 13 Call CalcLL recursive(LL, index + 1, current path, F path, Genetree)

Figure 12 :

 12 Figure 12: Posterior inference of female effective population size from 35 mtDNA samples from Yoruban individuals in the 1000 Genomes Project using our BESTT (first plot) from only 240 segregating sites and the BEAST BSP (second plot) from all the 15409 sites. Posterior median curves are depicted as solid black lines and 95% credible intervals by shaded regions.

%&'#( )*(+,(-.'

	A Data (Y)
		! " " # $ % % & & & '
	(	)	( * * ( * * * * * *
	)	)	( * * * ( * * * * *
	+	)	( * * * * ( ( * * *
	,	(	( * * * * * * ( ( (
	-	)	* * * * * * * * * *
	.	)	* * * * * * * * * *
	/	)	* ( ( * * * * * * *
	0	)	* * * * * * * * * *
	1	(	* * * * * * * * * *
		(.
	C Tajima's genealogy (g T )
		2				t 2
				3		t 3
	7	5	4	6		t 7 t 6 t 5 t 4
		8				t 8
	9					t 9
	12	16 10	13	15 11	14	t 16 t 15 t 13 t 12 t 14 t 11 t 10

Table 1 :

 1 Empirical measures of performance in the simulations described in the text

	Simulation	% ENV SRE MRW
	Instantaneous growth(n=10,s=31)	96	5.87 124352
	Instantaneous growth (n=10,s=63)	100	2.15 2296
	Instantaneous growth (n=10,s=120) 98	0.53 80
	Instantaneous growth (n=25)	90	0.40 3.43
	Instantaneous growth (n=35)	92	0.31 3.16
	Constant	100	0.30 1.16
	Exponential	100	0.35 5.45
	Exp. & const. (n=10, 1 locus)	100	4.31 22608
	Exp. & const. (n=10, 5 loci)	100	2.37 309.1
	Exp. & const. (n=10, 10 loci)	100	0.16 4.19

Table 2 :

 2 Performance comparison between BESTT and BEAST in simulations

			% ENV			SRE			MRW
	Simulation	BESTT EBSP Skygrid BESTT EBSP Skygrid BESTT EBSP Skygrid
	Constant	100	100	100	0.3	0.24	0.24	1.16	1.49	4.41
	Exponential 100	97	100	0.35	0.26	0.29	5.45	2.56	46.13
	Inst. growth 97	94	100	0.61	2.65	22.6	105.6	14.95 >1000

  BESTT and BEAST Comparison. Posterior inference from simulated data of n = 10 sequences under constant, exponential and instantaneous growth trajectories (rows) obtained from our method BESTT (first column), BEAST EBSP (second column) and BEAST Skygrid (third column). Shaded areas correspond to 95% credible intervals, solid black lines show posterior medians and dashed lines show the simulated truth. In the fourth column, we show the posterior density of the time to the most recent common ancestor (TMRCA) from the three methods: BESTT (black), BEAST EBSP (blue) and BEAST Skygrid (red). The true value of the TMRCA is depicted as a vertical dashed line. and Y * is the labeled data, in order to estimate (N (t)) t>0 . These two posterior distributions are the same when every individual of the sample has its own private mutation group and no shared mutation groups. Otherwise, the number of Tajima's trees compatible with observed data Y, i.e. Tajima's trees g T such that Pr(g T | Y) > 0, is smaller than the number of Kingman's trees compatible with observed labeled data Y

			Constant BESTT	Constant EBSP		Constant Skygrid			Constant TMRCA
																	-	truth
	N(t)	0.10 1.00			0.10 1.00						0.10 1.00			Density			---	BESTT EBSP Skygrid
		0.01			0.01						0.01
		0.7	0.5	0.3	0.1	0.6		0.4	0.2	0.0	0.6	0.4	0.2	0.0	0.3	0.5	0.7	0.9
			Exponential BESTT	Exponential EBSP		Exponential Skygrid		Exponential TMRCA
	N(t)	1e-02 1e+00 1e+02			1e-02 1e+00 1e+02						1e-02 1e+00 1e+02			Density
			0.8 0.6 0.4 0.2 0.0	0.8	0.6	0.4	0.2	0.0	0.8 0.6 0.4 0.2 0.0	0.5	0.7	0.9	1.1
		Instant growth BESTT	Instant growth EBSP		Instant growth Skygrid	Instant growth TMRCA
	N(t)	1e-03 1e+01			1e-02 1e+00 1e+02					N(t)	1e-03 1e+01
		0.08 0.06 0.04 0.02 0.00	0.08 0.06 0.04 0.02 0.00	0.08 0.06 0.04 0.02 0.00	0.04	0.06	0.08	0.10
				Time				Time				Time

Density

Figure 8: * (

Acknowledgements

The authors thank the editor and two anonymous referees whose suggestions considerably improved the manuscript. This research is supported in part by a National Institutes of Health grant R01-GM-131404 and the Alfred P. Sloan Foundation to J.A.P.. We want to acknowledge the developers of R-ape, R-phangorn and R-phylodyn that facilitated our implementations. A.V. was supported in part by the chaire program Modélisation Mathématique et Biodiversité of Veolia Environnement -École Polytechnique -Museum National d'Histoire Naturelle -Fondation X. A.V. and J.A.P. was supported by the France-Stanford Center for interdisciplinary Studies. This was work also

interval; that is, we choose i ∼ U ({n, n -1, . . . , 3}) and t * i ∼ U (t i+1 , t i-1 ). This is a symmetric proposal and the corresponding Metropolis-Hastings acceptance probability is

While these updates may seem unnecessary in light of the Hamiltonian updates of coalescence times (Section 2.8.3), we observed better performance of our MCMC sampler by including this additional proposal. One reason may be our choice of M in section 2.8.3 that does not account for the geometric structure of the posterior distribution of coalescent times. However, a better choice of M comes with higher computational burden than a simple local update of coalescent times.

Multiple Independent loci

Thus far, we have assumed our data consist of a single linked locus of s segregating sites. We can extend our methodology to l independent loci with s i segregating sites for i = 1, . . . , l. In this case, our data Y = (Y 1 , . . . , Y l ) consist of l aligned sequences with elements {0, 1}, where 0 indicates the ancestral allele as before. We then jointly estimate the Tajima's genealogies {g T i } l i=1 , precision parameter τ , and vector of log effective population sizes γ through their posterior distribution:

In Equation ( 14), we enforce that all loci follow the same effective population size trajectory but every locus can have its own mutation rate µ i .

Results

The performance of BESTT in applications to simulated data

We tested our new method, BESTT, on simulated data under four different demographic scenarios.

Note that in this section, N (t) is rescaled to the coalescent time scale, meaning that 1/N (t) is the pairwise rate of coalescence at time t in the past relative to the rate at the present time zero. We simulated genealogies under four different population size trajectories:

1. A period of exponential growth followed by constant size:

2. A trajector with instantaneous growth:

3. An exponential growth:

In order to leverage the constraints imposed by the data and the infinite-sites mutation model, we apply Dan Gusfield's perfect phylogeny algorithm [START_REF] Gusfield | Efficient algorithms for inferring evolutionary trees[END_REF] to represent sequence alignments as a gene tree. We exploit the gene tree representation for conditional likelihood calculations and for exploring the state space of ranked tree shapes.

For the calculation of the likelihood of the data conditioned on a given Tajima's genealogy, we augment the gene tree representation of the data with the Tajima's genealogy and map observed mutations to branches. We define a directed acyclic graph (DAG) with the augmented gene tree.

This new representation as a DAG allows for calculating the likelihood as a backtracking algorithm that transverses the gene tree from the leaves to the root. Our implementation's computational bottleneck lies in the likelihood calculation. Given a Tajima's genealogy, our likelihood algorithm sums over all possible allocation of mutation groups to branches. Although this number is generally much smaller than the number of labeled genealogies, our algorithm can be further optimized. In future studies, we will explore as sum-product type of algorithm for the likelihood calculation. In the present implementation we are able to infer effective population size trajectories from samples of size n ≈ 35 in a regular personal laptop computer within few hours.

Our statistical framework draws on Bayesian nonparametrics. We place a flexible geometric random walk process prior on the effective population size that allows us to recover population size trajectories with abrupt changes in simulations. The inference procedure proposed in this manuscript relies on Markov chain Monte Carlo (MCMC) methods with 3 large Gibbs block updates of: coalescent times, effective population size trajectory and ranked tree shape topology. We use Hamiltonian Monte Carlo updates for continuous random variables: coalescent times and effective population size; and a Metropolis Hastings sampler for exploring the space of ranked tree shapes.

For exploring the genealogical space, [START_REF] Markovtsova | The age of a unique event polymorphism[END_REF] suggest a joint local proposal for both coalescent times and topology. Here we restrict our attention to the topology alone. A future line of research includes the development of a joint local proposal of coalescent times and ranked tree shapes. We also envision that a joint sampler of coalescent times and effective population size trajectories should improve mixing and convergence.

Our method does not model recombination, population structure or selection. It assumes completely linked and neutral segments from individuals from a single population, and the infinite sites mutation model. While this model is a good approximation for some human molecular data, it is not appropriate for modeling molecular data from other organisms such as pathogens and viral populations. Finally, haplotype data of many organisms is usually sparse with few unique haplotypes presented at high frequencies. Since our algorithm exploits molecular data at the haplotype level, our proposed method is ideally suited for this scenario where the space of ranked tree shapes is drastically smaller than the space of labeled topolgies.

Appendix A Matrix representation of a ranked tree shape

Our algorithms exploit the following encoding of a ranked tree shape by a triangular matrix of size n × n, which we denote by F (Figure 10). Recall that, by convention, t n+1 = 0 and t 1 = +∞.

First, we declare that F i,j = 0 if j > i. Next, the number of lineages through time is encoded on the diagonal of F: F i,i = i for i in {2, 3, . . . , n}. Finally, for j < i, the entry F i,j denotes the number of lineages that do not coalesce in the time interval (t i+1 , t j ); in particular, F i,1 = 0 and for every i in {2, 3, . . . , n}, F n,i denotes the number of singletons (i.e., external branches that have not coalesced) in the time interval (t i+1 , t i ) (Figure 10). Other statistics of the ranked tree shape can be expressed in terms of the corresponding matrix F. Among them, the number c of cherries is equal to the number of times that the number of singletons decreases by 2 between lines i and i -1, since such an event means that the coalescence separating these two epochs was that of two external branches. That is,

Appendix B

Detailed allocation of mutation groups along g T . The latent allocation random variables {A j } are constrained by the information in the Tajima's genealogy g T . In a given g T , every subtree is labeled by its ranking from past to present (Figure 10). Subtree i is subtended by branch b i with length l i , for i = 2, . . . , n. We will assume that l 2 , the length of the root branch, is 0. Let c be the number of cherries (nodes with two leaves) in g T ; the two branches of a given cherry share the same label b j ∈ {b n+1 , . . . , b n+c }. The actual label of external branches is arbitrary but, for ease of exposition in our figures, we first label the cherries' branches from left to right by {b n+1 , . . . , b n+c }; singleton branches are labeled from left to right by b n+c+1 , . . . , b 2n-c (Figure 2C). As mentioned before, the length of X j is the number of the corresponding sister nodes in T that were grouped together in forming node Z j . In this case, A j = (A j,1 , . . . , A j,|ch(j)| ) denotes a collection of |ch(j)| vectors of branch labels in g T subtending the child-node subtrees of node Z j . A j,1 corresponds to the branch subtending from the leftmost child node of Z j on D, A j,2 corresponds to the branch subtending from the next child node of Z j , etc., and A j,|ch(j)| corresponds to the branch subtending from the rightmost child node of Z j on D. Observe that, since we group some of the leaf nodes in T into a single node in D, any A j,k may be a vector of branch labels; for example A 1,1 = (b 12 , b 9 ) and A 1,2 = b 10 in Figure 3B.

Appendix C

Algorithms for conditional likelihood calculation. The following two algorithms detail the calculation of Pr(Y | g T , µ). Y is encoded in GeneT ree, the observed data as a Tree structure.

Each node in GeneT ree has number of descendants (or lineages) and mutation information attached to it. Tajima's genealogy g T is encoded as F path that contains the ranked tree shape F n and times that contains the vector of coalescent times t multiplied by the mutation rate µ.

Algorithm 3 Transition proposals for ranked tree shapes

Input:

2. Sample with uniform discrete probability a coalescent event k from the set {3, . . . , n}. Set

If the lineages coalescing at time t k are singletons (Figure 4, Case A, lineages 1 and 2 in F n )) (a) No sampling required to distinguish between two singletons. Set q 2 = 1.

(b) Update F * n : merge one singleton with the lineage coalescing at k -1 (excluding lineage k) in F n , then merge the second singleton at time t k-1 with lineage k.

(c) Compute the probability q 2 of restoring the ordering of F * n to F n .

Else

(a) Sample one the lineages coalescing at time t k with uniform discrete probability. Set q 2 = 1 2 (b) Update F * n : merge the sampled lineage with the one coalescing at time t k-1 in F n . At time t k-1 , merge the lineage not sampled with the new lineage k.

(c) Compute the probability q 2 of restoring the ordering of F * n to F n .

Else (Figure 4, Case B)

Swap the coalescent events. Lineages descending from k are now set to coalesce at time t k-1 and lineages previously descending from k -1 are now set to coalesce at time t k . Set q 2 = 1 and q 2 = 1.

Appendix D

We replicated the BEAST EBSP Analysis of the 35 Yoruban individuals from the 1000 Genomes Project phase 3 using the whole mtDNA coding region consisting of 15409 sites. In both cases we assumed the Jukes-Cantor mutation model [START_REF] Jukes | Evolution of protein molecules[END_REF]. Figure 11 shows the comparison between EBSP inference from the 240 segregating sites retained in section 4 that are compatible with the infinite sites mutation model assumption. In both cases we recover very similar trajectories.

In addition, we compared our results with BEAST Bayesian Skyline Plot (BSP) [START_REF] Drummond | Reconstructing genealogies of serial samples under the assumption of a molecular clock using serial-sample UPGMA[END_REF]. For our reduced dataset of 240 segregating sites, we could not generate valid inference of N (t) with Metropolis-Hastings acceptance probability greater than 0. Instead we were able to generate results with BEAST BSP from the complete dataset of 15409 sites. The comparison of our method from 240 segregating sites to BEAST BSP from 15409 sites is depicted in Figure 12.

Appendix E

In Figure 13A, we show the data from Figure 2A with an additional haplotype (10) with frequency 1 and an additional column grouped with mutation group h (not shown in the table). In 13B we show the corresponding perfect phylogeny. This new perfect phylogeny has a new tip with black label 1 (frequency) subtending from a branch with 0 mutations (red label). The path from the leaf