Genomic data analysis in tree spaces

In this subsection we quickly explain the foundations of metric geometry. See [7, 8] for comprehensive textbooks on the subject. A metric space $(X,d)$ is a set $X$ equipped with a distance function $d\colon X\times X\to\mathbb{R}^{\geq 0}$ having the properties that $d(x,y)=d(y,x)$, $d(x,y)=0$ if and only if $x=y$, and $d(x,z)\leq d(x,y)+d(y,z)$ for all $x,y,z\in X$ (the triangle inequality). Although metric spaces often arise in contexts in which there is not an evident notion of geometric structure, it turns out that under very mild hypotheses a metric space $(X,d)$ can be endowed with structures analogous to those arising on Riemannian manifolds. A metric space is a length space if the distance $d(x,y)$ is realized as the infimum of the lengths of paths joining $x$ and $y$. A length space $(X,d)$ is a geodesic metric space if any two points $x$ and $y$ can be joined by a path with length precisely $d(x,y)$. A key insight of Alexandrov is that curvature makes sense in any geodesic metric space [1].

The idea is that the curvature of a space can be detected by considering the behavior of the area of triangles, and triangles can be defined in any geodesic metric spaces. Specifically, given points $p,q,r$, we have the triangle $T=[p,q,r]$ with edges the paths that realize the distances $d(p,q)$, $d(p,r)$, and $d(q,r)$. The connection between curvature and area of triangles comes from the observation that given side lengths $(\ell_{1},\ell_{2},\ell_{3})\subset\mathbb{R}^{3}$, a triangle with these side lengths on the surface of the Earth is “fatter” than the corresponding triangle on a Euclidean plane. To be precise, we consider the distance from a vertex of the triangle to a point $p$ on the opposite side — in a fat triangle, this distance will be larger than in the the corresponding Euclidean triangle (and smaller in a thin triangle).

Given a triangle $T=[p,q,r]$ in $(X,d)$, we can find a corresponding triangle $\tilde{T}$ in Euclidean space with the same edge lengths. Given a point $z$ on the edge $[p,q]$, a comparison point in $\tilde{T}$ is a point $\tilde{z}$ on the corresponding edge $[\tilde{p},\tilde{q}]$ such that $d_{E}(\tilde{z},\tilde{p})=d(z,p)$ and $d_{E}(\tilde{z},\tilde{q})=d(z,q)$. (Where here $d_{E}$ denotes the Euclidean metric.) We say that a triangle $T$ in $M$ satisfies the $\textrm{CAT}(0)$ inequality if for every such pair $(z,\tilde{z})$, we have $d(r,z)\leq d_{E}(\tilde{r},\tilde{z})$. If every triangle in $M$ satisfies the $\textrm{CAT}(0)$ inequality then we say that $M$ is a $\textrm{CAT}(0)$ space.

More generally, let $M_{\kappa}$ denote the unique two-dimensional Riemannian manifold with curvature $\kappa$. The diameter of $M_{\kappa}$ will be denoted $D_{\kappa}$. Then we say that a geodesic metric space $M$ is $\textrm{CAT}(\kappa)$ if every triangle in $M$ with perimeter $\leq 2D_{\kappa}$ satisfies the inequality above for the corresponding comparison triangle in $M_{\kappa}$. If $\kappa^{\prime}\leq\kappa$, any $CAT(\kappa^{\prime})$ space is also $\textrm{CAT}(\kappa)$. A $n$-dimensional Riemannian manifold $M$ that is sufficiently smooth has sectional curvature $\leq\kappa$ if and only if $M$ (regarded as a metric space) is $\textrm{CAT}(\kappa)$. For example, Euclidean spaces are $\textrm{CAT}(0)$, spheres are $\textrm{CAT}(1)$, and hyperbolic spaces are $\textrm{CAT}(-1)$.

As described, $\textrm{CAT}(\kappa)$ is a global condition; we will say that a metric space $(X,d)$ is locally $\textrm{CAT}(\kappa)$ if for every $x$ there exists a radius $r_{x}$ such that $B_{r_{x}}(x)\subseteq X$ is $\textrm{CAT}(\kappa)$. For example, the flat torus (obtained by identifying opposite edges in a rectangle) is locally $\textrm{CAT}(0)$ but not globally $\textrm{CAT}(0)$. The Cartan-Hadamard theorem implies that a simply-connected metric space that is locally $\textrm{CAT}(0)$ is also globally $\textrm{CAT}(0)$.

A remarkably productive observation of Gromov is that many geometric properties of Riemannian manifolds are shared by $\textrm{CAT}(\kappa)$ spaces. In particular, $\textrm{CAT}(\kappa)$ spaces with $\kappa\leq 0$ (referred to as non-positively curved metric spaces) admit unique geodesics joining each pair of points $x$ and $y$, balls $B_{\epsilon}(x)$ are convex and contractible for all $x$ and $\epsilon\geq 0$, and midpoints of geodesics are well-behaved. As a consequence, there exist well-defined notions of mean and variance of a set of points, and more generally one can develop some of the foundations of classical statistics, as we review below in Section 3.

In this subsection, we review the theory of cubical complexes, which provide a rich source of examples of $\textrm{CAT}(0)$ metric spaces (again, see [8] or [7] for textbook treatments). It is in general very difficult to determine for an arbitrary metric space whether it is $\textrm{CAT}(\kappa)$ for any given $\kappa$. Even for finite polyhedra where the metric is induced from the Euclidean metric on each face, this problem does not have a general solution. The important of cubical complexes in this context comes from an effective criterion for determining if they are non-positively curved (i.e., $\textrm{CAT}(0)$).

Let $I^{n}\subseteq\mathbb{R}^{n}$ denote the $n$-dimensional unit cube $[0,1]\times\ldots\times[0,1]$, regarded as inheriting a metric structure from the standard metric on $\mathbb{R}^{n}$. A codimension $k$ face of the cube $I^{n}$ is determined by fixing $k$ coordinates to be in the set $\{0,1\}$. A cubical complex is a metric space obtained by gluing together cubes via the data of isometries of faces, subject to the condition that two cubes are connected by at most a single face identification and no cube is glued to itself. The metric structure is the length metric induced from the Euclidean metric on the cubes, i.e., the distance between $x$ and $y$ is the infimum of the lengths over all paths from $x$ to $y$ that can be expressed as the union of finitely many segments each contained within a cube. When the cubical complex $C$ is finite or locally finite, results of Bridson [6] and Moussong [33] imply that $C$ is a complete geodesic metric space.

Gromov gave a criterion for a cubical complex to be $\textrm{CAT}(0)$ that is often possible to check in practice. In order to explain this criterion, we need to review the notion of the link of a vertex in a cubical complex.

Fix a vertex $v$ in a cubical complex $C$ and a cube $C_{i}\cong I^{m}\subseteq C$ such that $v$ is a vertex of $C_{i}$. For fixed $\epsilon>0$, the all-right spherical simplex associated to $(C_{i},v)$ is the subset

The set $S(C_{i},v)$ has a metric induced by the Euclidean angle metric. The faces of $S(C_{i},v)$ are defined as the intersections of $S(C_{i},v)$ with faces of $C_{i}$; equivalently, these are the all-right spherical simplexes associated to faces of $C_{i}$. The collection of all-right spherical simplices for all pairs $(C_{i},v)$ forms a polyhedral complex with metric given by the length metric induced from the angle metrics; this is referred to as a spherical complex. Forgetting the metric structure, the all-right spherical simplices also form an abstract simplicial complex. (Recall that an abstract simplicial complex is simply a set of subsets of a set $V$ that is closed under passage to subsets.)

The link of a vertex $v$ in a cubical complex $C$ is the spherical complex obtained as the subset

for fixed $0<\epsilon<1$. Gromov’s criterion now states that the cubical complex $C$ is locally $\textrm{CAT}(0)$ if and only if the link is $\textrm{CAT}(1)$ or the abstract simplicial complex underlying the link is flag. (Recall that a flag complex is a simplicial complex in which a $k$-simplex is in the complex if and only its $1$-dimensional faces are in the complex.)

As an easy application of Gromov’s criterion, we conclude the section by showing the standard result that the Cartesian product of locally $\textrm{CAT}(0)$ cubical complexes is itself a locally $\textrm{CAT}(0)$ cubical complex. Let $X$ and $Y$ be cubical complexes that are $\textrm{CAT}(0)$. Since $I^{n}\times I^{m}\cong I^{n+m}$, it is clear that $X\times Y$ has the structure of a cubical complex. The set of vertices of $X\times Y$ is given by the product of the sets of vertices of $X$ and $Y$ respectively. The link of a vertex $(v,v^{\prime})$ in $X\times Y$, regarded as an abstract simplicial complex, is the join of $L(v)$ and $L(v^{\prime})$, which we denote $L(v)\ast L(v^{\prime})$ (the join of complexes $S_{1}$ and $S_{2}$ is obtained by considering all pairwise unions of elements of $S_{1}$ and $S_{2}$). Finally, since the join of flag complexes is easily seen to be a flag complex, Gromov’s criterion now implies that $X\times Y$ is locally $\textrm{CAT}(0)$.

A phylogenetic tree with $m$ leaves is a weighted, connected graph with no cycles, having $m$ distinguished vertices (referred to as leaves) of degree $1$ and labeled $\{1,\ldots,m\}$. All the other vertices are of degree $\geq 3$. We refer to edges that terminate in leaves as external edges and the remaining edges are internal.

The space $\textrm{BHV}_{m}$ of isometry classes of rooted phylogenetic trees with $m$-labelled leaves where the nonzero weights are on the internal branches was introduced and studied by Billera, Holmes, and Vogtmann [5]. The space $\textrm{BHV}_{m}$ is constructed by gluing together $(2m-3)!!$ positive orthants $\mathbb{R}^{m}_{\geq 0}$; each orthant corresponds to a particular tree topology, with the coordinates specifying the lengths of the edges. A point in the interior of an orthant represents a binary tree; if any of the coordinates are $0$, the tree is obtained from a binary tree by collapsing some of the edges. We glue orthants together such that a (non-binary) tree is on the boundary between two orthants when it can be obtained by collapsing edges from either tree geometry. Put another way, two tree topologies are adjacent when they are connected by a rotation, i.e., one topology can be generated from the other by collapsing an edge to length $0$ and then expanding out another edge from the incident vertex.

The metric on $\textrm{BHV}_{m}$ is induced from the standard Euclidean distance on each of the orthants. For two trees $t_{1}$ and $t_{2}$ which are both in a given orthant, the distance $d_{\textrm{BHV}_{m}}(t_{1},t_{2})$ is defined to be the Euclidean distance between the points specified by the weights on the edges. For two trees which are in different quadrants, there exist (many) paths connecting them which consist of a finite number of straight lines in each quadrant. The length of such a path is the sum of the lengths of these lines, and the distance $d_{\textrm{BHV}_{m}}(t_{1},t_{2})$ is then the minimum length over all such paths. For many points, the shortest path goes through the “cone point”, the star tree in which all internal edges are zero.

Allowing potentially nonzero weights for the $m$ external leaves corresponds to taking the cartesian product with an $m$-dimensional orthant. We will focus on the space

which we refer as the evolutionary moduli space (the $m-1$ index arises from the fact that we consider unrooted trees.) There is a metric on $\Sigma_{m}$ induced from the metric on $\textrm{BHV}_{m-1}$. Specifically, for a tree $t$, let $t(i)$ denote the length of the external edge associated to the vertex $i$. Then

where $\bar{t}_{i}$ denotes the tree in $\textrm{BHV}_{m-1}$ obtained by forgetting the lengths of the external edges (e.g., see [39]). As explained in  [5, §4.2], efficiently computing the metric on $\Sigma_{m}$ is a nontrivial problem, although there exists a polynomial-time algorithm [39].

The main result of Billera, Holmes, and Vogtmann is that the length metric on $\textrm{BHV}_{n}$ endows this space with a (global) $\textrm{CAT}(0)$ structure. By subdividing each orthant into cubes in the evident fashion, $\Sigma_{m}$ is naturally a cubical complex where the metric we have described is the one induced from the Euclidean metric on the cubes; a straightforward combinatorial analysis of the link of $\Sigma_{m}$ implies the result via Gromov’s criterion. In addition, $\Sigma_{m}$ is clearly a complete and separable metric space; any tree can be approximated by a sequence of trees in the same orthant that have rational edge lengths.

In evolutionary applications, we are often interested in classifying and comparing distinct behaviors by understanding the relative lengths of edges: rescaling edge lengths should not change the relationship between the branches [57]. Motivated by this consideration, we define $\mathbb{P}\Sigma_{m}$ to be the subspace of $\Sigma_{m}$ consisting of the points $\{t_{i}\}$ in each orthant for which the constraint $\sum_{i}t_{i}=1$ holds.

We denote the space of trees with internal edges of fixed length by $\tau_{m-1}$. The space of $m$ external branches whose lengths sum to $1$ is the standard $m-1$ dimensional simplex $\Delta_{m-1}$ in $\mathbb{R}^{m}$. The constraint that the length of internal branches plus the external branches sum to $1$ implies that

where here $\star$ denotes the join of two spaces, using the Milnor model of the join. We can also describe $\mathbb{P}\Sigma_{m}$ as the link on the origin in $\Sigma_{m}$.

There are various possible natural metrics to consider on $\mathbb{P}\Sigma_{m}$. The simplest way to endow $\mathbb{P}\Sigma_{m}$ with a metric is to use the induced intrinsic metric specified by paths in $\Sigma_{m}$ constrained to lie entirely within $\mathbb{P}\Sigma_{m}$. From the perspective of metric geometry, the characterization of $\mathbb{P}\Sigma_{m}$ as the link of the origin endows it with a “spherical” metric, and Gromov’s criteria imply that with this metric, $\mathbb{P}\Sigma_{m}$ is a $\textrm{CAT}(1)$ space. (Alternatively, $\mathbb{P}\Sigma_{m}$ is the spherical join of $\tau_{m-1}$ and the spherical realization of the $\Delta_{m}$; since $\tau_{m-1}$ and $\Delta_{m}$ are $\textrm{CAT}(1)$, so is their spherical join [7, II.3.15].) The theory of polyhedral complexes implies that in either case $\mathbb{P}\Sigma_{m}$ is a complete geodesic metric space [7, I.7.19], and it is evidently separable.

Moreover, with the induced intrinsic metric $\mathbb{P}\Sigma_{m}$ is in fact a $\textrm{CAT}(0)$ space; although $\tau_{m-1}$ has points which are not connected by unique geodesics (see Section 2.5 below for a more detailed discussion), the join with $\Delta_{m}$ introduces a new “cone direction” that changes the geometry.

To compute the intrinsic metric on $\mathbb{P}\Sigma_{m}$, we use $\epsilon$-nets and a local-to-global construction. Recall that a set of points $S$ in a metric space $(X,\partial)$ is an $\epsilon$-net if for every $z\in X$, there exists $q\in S$ such that $\partial(z,x)<\epsilon$. For a compact metric space equipped with a probability measure such that all non-empty balls in the metric space have nonzero measure, we can produce an $\epsilon$-net by sampling. More precisely, it is straightforward to show that given a finite collection of measurable sets $\{A_{1},A_{2},\ldots,A_{k}\}$ and a probability measure $\mu$ on $\cup_{i}A_{i}$ such that $\mu(A_{i})\geq\alpha>0$, then given at least

samples, with probability $1-\delta$ there is at least one sample in every $A_{i}$ [35, 5.1].

Next, suppose that we have a metric space $(X,\partial)$ where there exists a constant $\kappa$ such that if $\partial(x,y)<\kappa$, it is easy to compute $\partial(x,y)$. An algorithm for approximating $\partial$ on all of $X$ is then to take a dense sample $S\subset X$, form the graph $G$ with vertices the points of $S$ and edges between $x$ and $y$ when $\partial(x,y)<\kappa$, and define the distance between $x$ and $y$ in $X$ to be the graph metric on $G$ between the nearest points to $x$ and $y$ in $S$. This distance can be efficiently computed using Dijkstra’s algorithm [12].

When $S$ is an $\epsilon$-net for sufficiently small $\epsilon$ relative to $\kappa$, we can describe the quality of the resulting approximation to $\partial$ [4, Thm. 2]. In particular, if $\epsilon<\frac{\kappa}{4}$, then

Putting this all together, to approximate the metric on $\mathbb{P}\Sigma_{m}$ we take the union of $\epsilon$-nets on all of the simplices (including the faces) and form a $\kappa$-approximation $\partial_{G}$ as above. In practice, the required density of samples is determined by looking at when the approximation converges (i.e., when the change in distances drops below a specified precision bound). We sample densely on each simplex representing a tree topology on $m$ leaves, and explicitly include certain key singular points of the projective space. Figure 3 contains two visualizations of the projective space, using force-directed layouts of the $\epsilon$-nets constructed on 4–leaved and 5–leaved trees respectively.

We now describe the size of $\mathbb{P}\Sigma_{m}$ (see also [5, 3.3] for a related discussion). For simplicity, we will focus on the link in $\textrm{BHV}_{n}$, which we will denote by $\mathcal{P}_{n}$, and temporarily ignore the join with the simplex coming from the external edge lengths. Observe that adjacent top-dimensional simplices in $\mathcal{P}_{n}$ differ by a rotation of tree topologies, where a rotation collapses an internal edge and then expands out from the resulting node. Next, recall that the homotopy type of $\mathcal{P}_{n}$ is a wedge of $(n-1)!$ spheres of dimension $(n-3)$ [44] (and see also [11, Thm. 6]). Moreover, we can explicitly describe these spheres, as follows.

As discussed in [11, Prop. 1] and [5, §3.1], the boundary of the dual polytope to standard associahedron on $n$ letters (parametrizing parenthesizations of $n$ terms) embeds in many different ways into $\mathcal{P}_{n}$. Following [11], let us denote this boundary by $\mathcal{K}_{n}$. Explicitly $\mathcal{K}_{n}$ is a simplicial sphere of dimension $(n-3)$ where a $k$-simplex corresponds to a planar rooted tree with $n$ leaves and $k+1$ internal edges. Then the homotopy type of $\mathcal{P}_{n}$ can be described in terms of various embedded copies of $\mathcal{K}_{n}$. As a consequence, to understand the size of $\mathcal{P}_{n}$, we need to compute the diameter of $\mathcal{K}_{n}$.

For convenience, we describe this diameter in terms of counts of simplices; the actual value can then be obtained by multiplying by the diameter of a simplex. In this guise, the problem is an old one which can be described in many different forms, perhaps most relevantly as the computation of maximal rotation distances between binary trees.

The main result here is that, for unrooted trees on $n$ leaves, the diameter of $\mathcal{K}_{n}$ is $2n-8$ for $n>11$ [41]; this bound was established asymptotically (for sufficiently large but indeterminate $n$) in [49]. For smaller values, we have the following table (taken from [49, §2.3]) of explicit values:

More generally, the maximum rotation distance between labelled trees on $n$ leaves is $O(n\log n)$ [50]. Of course, the cone point associated to the join with standard simplex means that the size is considerably smaller.

In order to study probability distributions in evolutionary moduli spaces, it is necessary to have reasonable notions of moments of the distribution, expectation of random variables, and analogues of the law of large numbers. Since $\Sigma_{m}$ and $\mathbb{P}\Sigma_{m}$ are $\textrm{CAT}(0)$ spaces, points are connected by unique geodesics and there is a sensible notion of a centroid of a collection of points. Discussion of statistical inference in $\Sigma_{m}$ was initiated in [5], and subsequently Holmes has written extensively on this topic [27, 27, 28] (and see also [18]). More generally, Sturm explains how to study probability measures on general $\textrm{CAT}(0)$ spaces [53]. He shows that there are reasonable notions of moments of distribution, expectation of random variables, and analogues of the law of large numbers on $\textrm{CAT}(0)$ spaces.

Sturm provides an iterative procedure for computing the mean and variance of a set of points in $\Sigma_{m}$, and by exploiting the local geometric structure of $\Sigma_{m}$, Miller, Owen, and Provan produce somewhat more efficient algorithms for computing the mean [32]. Furthermore, Sturm proves versions of Jensen’s inequality and the law of large numbers in this context. The situation for the central limit theorem is less satisfactory. Barden, Le, and Owen study central limit theorems for Fréchet means in $\Sigma_{m}$ [3]; as they explain, the situation exhibits non-classical behavior and the limiting distributions depend on the codimension of the simplex in which the mean lies. Finally, there has been some work on principal components analysis (PCA) in $\Sigma_{m}$ [37].

However, in contrast to classical statistics on $\mathbb{R}^{n}$, we do not know many sensible analytically-defined distributions on the evolutionary moduli spaces. Billera-Vogtmann-Holmes briefly introduce a family of Mallows distribution on $\Sigma_{m}$ with density function

for fixed $t_{1}\in\Sigma_{m}$, and an analogous family can be defined on $\mathbb{P}\Sigma_{m}$. Sampling from these distributions is not easy; in general, the behavior of distributions on $\Sigma_{m}$ and sampling algorithms is somewhat perverse due to the pathological behavior near the origin due to the exponential growth in the mass of an $\epsilon$ ball. A much more tractable source of distributions on $\Sigma_{m}$ and $\mathbb{P}\Sigma_{m}$ arise from resampling from a given set of empirical data points.

One way to grapple with the difficulties in dealing with distributions on $\Sigma_{m}$ and $\mathbb{P}\Sigma_{m}$ is to instead study associated projections into distributions on Euclidean space. The advantage of this approach is evident; we are now in a setting where the theory of moments, the central limit theorem, and asymptotic consistency for resampling procedures are all very familiar. Of course, it is important to keep in mind that the moments derived in this setting will reflect the geometry of the evolutionary moduli space is complicated ways, and inverses to the projections will not usually exist. Nonetheless, for purposes of many kinds of statistical tests (e.g., hypothesis testing about distributions generating observed samples), this approach can be very effective. There are a number of natural ways to map metric measure spaces into Euclidean space; in this section, we discuss several strategies derived from the use of the metric.

An intrinsic map comes from looking at the distance distribution on $\mathbb{R}$ induced by $\partial_{M}$. Specifically, given a Borel distribution $\Psi$ on $(M,\partial_{M})$, the product distribution $\Psi\otimes\Psi$ on $M\times M$ induces a distibution on $\mathbb{R}$ via $\partial_{M}$. Applying this construction to the empirical measure on a finite sample yields the empirical distance distribution. More generally, for any fixed $n$, we can consider the distribution on $\mathbb{R}^{n^{2}}$ induced by taking the product distribution $\Psi^{\otimes n}$ on $M^{\times n}$ and applying $\partial_{M}$ to produce the $n\times n$ matrix of distances. Gromov’s “mm-reconstruction theorem” showed that in the limit as $n\to\infty$, the distance matrix distributions completely characterize the distribution $\Psi$ on $(M,\partial_{M})$ [23]. Once again, given a sufficiently large finite sample, we can construct the empirical distance matrix distributions for any fixed $n$.

Another approach involves choosing a fixed set of $n$ landmarks and considering the vector of distances from a fixed point to the landmarks. Given a set $L=\{t_{1},t_{2},\ldots,t_{n}\}\subset M$, there is a continuous map

specified by

Pushforward along $d_{L}$ again induces a distribution on $\mathbb{R}^{n}$ from one on $M$. One expects that as $k$ increases (provided the landmarks are “generic”), the induced distributions in $\mathbb{R}^{n}$ will characterize the distribution on $M$.

In both cases, choice of the parameter $k$ depends on some sense of the intrinsic dimension of the support of the distribution as well as the number of points available (in the case of finite samples). Unfortunately, the required $k$ may well be quite large.

Finally, there is a substantial body of work on low-distortion embeddings of finite metric spaces into $\ell^{p}$ spaces, in particular Euclidean spaces. Recall that the distortion of a non-contractive (distance expanding) embedding of metric spaces $f\colon(X,\partial_{X})\to(Y,\partial_{Y})$ is given by $\sup_{x_{1}\neq x_{2}}\frac{\partial_{Y}(f(x_{1}),f(x_{2}))}{\partial_{X}(x_{1},x_{2})}$. Notably, Abraham, Bartal, and Neiman show that one can construct a probabilistic embedding of an $n$-point finite metric space into an $O(\log n)$ dimensional space with distortion $O(\log n)$. Pushing forward distributions on $\Sigma_{m}$ and $\mathbb{P}\Sigma_{m}$ provide another way of reducing statistical questions to Euclidean space.

Given a set of samples $X\subset\Sigma_{m}$ and a partition $X=\mathcal{C}_{1}\cup\mathcal{C}_{2}\ldots\cup\mathcal{C}_{n}$ (where $\mathcal{C}_{i}\cap\mathcal{C}_{j}=\emptyset$), it is often useful to be able to determine whether or not the different $\mathcal{C}_{i}$ were generated from the same or different underlying distributions. For instance, $\mathcal{C}_{1}$ might represents samples from patients who received treatment and $\mathcal{C}_{2}$ is untreated patients, or the different groups $\mathcal{C}_{i}$ represent different observed genetic markers. Based on the discussion of the previous two subsections, we can study this problem directly in $\Sigma_{m}$ or in via projections to $\mathbb{R}^{n}$.

Te Fréchet mean and variance provides a summary of each collection of samples $\mathcal{C}_{i}$. A standard comparison between groups is then given by the distance

between the means. In order to understand the variability due to sampling, we can use bootstrap resampling (or more general $k$ out of $n$ resampling without replacement) to generate confidence intervals for the value of $\theta_{i}$. Asymptotic consistency for the bootstrap follows from the fact that the VC dimension of the collections of balls in $\Sigma_{m}$ and $\mathbb{P}\Sigma_{m}$ is bounded, via the usual criteria [20, 21, 22].

However, it is often simpler to consider tests induced by the projections into $\mathbb{R}^{n}$ discussed above. Here, we can compare collections $\mathcal{C}_{i}$ and $\mathcal{C}_{j}$ by using any of the many standard non-parametric comparison techniques for real distributions, for example $\chi^{2}$ tests or two-sample Kolmogorov-Smirnov tests.

One pervasive problem in clinical applications is that often the number of samples is quite small, and so we are often far from the asymptotic regime for statistical tests. Standard small-sample corrections can be applied. However, for this reason a machine learning approach to analyzing the data is often more useful.

Given the difficulties with statistical inference in $\Sigma_{m}$ and $\mathbb{P}\Sigma_{m}$, it is useful to complement these approaches with techniques from machine learning. The most basic family of techniques we might consider is clustering, a kind of unsupervised learning. Here, given a finite set $X$ in $\Sigma_{m}$ or $\mathbb{P}\Sigma_{m}$, we search for a partitionings of the points into clusters which optimize some criterion for the “goodness” of the clustering.

Regarding $\Sigma_{m}$ and $\mathbb{P}\Sigma_{m}$ simply as metric spaces, we can apply standard clustering algorithms that operate on arbitrary metric spaces. For example, we can apply standard $k$-means clustering, using the centroids as defined above. A related alternative is the $k$-medoids algorithm. Like $k$-means, $k$-medoids seeks partitions which are optimal in the sense of minimizing the sum of squared distances; the cost function for a cluster $C=\{x_{i},x_{j}\}$ is given by $\sum_{i<j}\partial(x_{i},x_{j})^{2}$. But instead of using cluster centroids as in $k$-means, cluster assignments are determined by medoids, which are points $z\in C$ that minimize $\sum_{i}d(z,i)$. The advantage of $k$-medoids over $k$-means is that the problem of finding a centroid in $\Sigma_{m}$ or $\mathbb{P}\Sigma_{m}$ is avoided.

Another natural family of clustering algorithms comes from spectral clustering techniques. Recall that spectral clustering can be applied to finite subsets of any metric space; one constructs an embedding into Euclidean space using the graph Laplacian associated to a graph encoding the local metric structure of the set of points and then performs $k$-means clustering. As such, spectral clustering can be applied both to $\Sigma_{m}$ and $\mathbb{P}\Sigma_{m}$. However, as illustrated in Figure 4, there is substantial distortion under such embeddings for low dimensions.

Although clustering algorithms are very useful for exploratory data analysis, for clinical applications we expect that classification problems are more salient. Specifically, a temporal sequence of tumor samples will be linked with a categorical or numeric label denoting the clinical management of the patient. We would then like to predict patient outcomes or expect response to treatment using a discriminative supervised learning algorithm operating in $\Sigma_{m}$ or $\mathbb{P}\Sigma_{m}$. Analogous to our use of $k$–medoids clustering for unsupervised grouping, the most basic algorithm for supervised learning is a $k$–nearest neighbor ($k$–NN) predictor. In this algorithm the predicted label for a given point is generated by taking a majority or weighted vote over the labels of the $k$ nearest trees. The optimal value of $k$ then specifies an order–$k$ Voronoi tesselation of the space that provides a description of the sizes of the predictive neighborhoods surrounding each element of the data set. We use this classification algorithm to study clinical correlates of trees determined by tumor samples from glioma patients; see Figure 9.

Let $\mathcal{E}_{m}$ denote either $\Sigma_{m}$ or $\mathbb{P}\Sigma_{m}$.

A representative example of $\Psi_{S}$ is shown in Figure 5.

Using $\Psi$, we can describe a number of dimensionality reduction procedures. The most basic example is simply to exhaustively subsample the labels. Let $\mathcal{D}(\Sigma_{m})$ denote the set of distributions on $\Sigma_{m}$.

(In practice, we approximate $\Psi_{k}$ using Monte Carlo approximations.)

Often there is additional structure in the labels that can be exploited. For instance, in many natural examples, the genomic data has a natural chronological ordering. When this holds, sliding windows over the labels induces an ordering on subtrees generated by $\Psi_{S}$. Rather than just regarding such a sequence as a distribution, the ordering makes it sensible to consider the associated trees as forming a piecewise-linear curve in $\Sigma_{k}$. (Note that given a set of points in $\Sigma_{k}$ it is always reasonable to form the associated piecewise-linear curve because each pair of points is connected by a unique geodesic.) Let $\mathcal{C}_{k}(\Sigma_{m})$ denote the set of piecewise linear curves in $\Sigma_{m}$; equivalently, $\mathcal{C}_{k}(\Sigma_{m})$ can be thought as the set of ordered sequences in $\Sigma_{m}$ of cardinality $k$. A schematic example of this sequential operation is given in Figure 6.

There are many variants of $\Psi_{C}$ depending on the precise strategy for windowing that is employed.

We have described tree dimensionality reduction in terms of the map $\Psi$, which is an operation on tree spaces. In practice, this technique would be applied by producing a very large phylogenetic tree from the raw data and subsequently applying $\Psi$. However, there is an alternative form of tree dimensionality reduction that instead subsamples the raw data to produce smaller phylogenetic trees. In this section, we discuss the relationship between these two procedures in the context of neighbor-joining.

Neighbor-joining is an algorithm for producing a tree from metric data (or more broadly, a dissimilarity measure) [17]. That is, the input is a set of points $X$ and a metric $\partial_{X}\colon X\times X\to\mathbb{R}$. One of the main theorems about consistency of neighbor-joining is that when $\partial_{X}$ is “close” to a tree metric $\partial_{T}$, neighbor-joining recovers $T$. Recall that given a tree $T$, the associated metric $\partial_{T}$ is defined by taking the distance between leaves $i$ and $j$ to be

where $P_{ij}$ is the unique path in $T$ from $i$ to $j$ and $\ell(e)$ is the length of the edge $e$.

The specific consistency theorem we use is due to Atteson [2]: if

then neighbor-joining recovers $T$. In this case we say that $(X,\partial_{X})$ is consistent with $T$.

For a metric space $(X,\partial_{X})$, let $T(X)$ denote the tree obtained from neighbor-joining applied to $(X,\partial_{X})$. As a consequence of Proposition 4.4, given a metric space $(X,\partial_{X})$ that is consistent with $T(X)$, the distribution $\Psi_{k}(T(X))$ is identical to the distribution $\{T(X_{S})\}$ where $S$ varies over all subsets of $X$ of cardinality $k$ and $X_{S}$ denotes the metric space structure on $S$ induced by $\partial_{X}$.

In order to apply tree dimensionality reduction in the face of potentially noisy data, we would like to know that small random perturbation of the original sample results in a distribution of subsamples that is “close” in some sense (e.g., small shifts in the centroid in $\Sigma_{m}$). Conversely, if two distributions of subsamples are suitably close, we would like to be able to conclude that the sampled trees are also close.

We begin with a lemma describing the interaction of $\Psi_{S}$ and the boundaries of orthants.

In light of Lemma 4.6, the projection $\Psi_{S}$ preserves paths. The other thing we need to understand is the potential increase in length caused by applying $\Psi_{S}$. Specifically, we need to consider the impact of the addition of edge lengths that occurs when a degree 2 vertex is produced by the reduction process. In the simplest case, we are considering the map $\mathbb{R}^{2}\to\mathbb{R}$ specified by $(x_{1},x_{2})\mapsto x_{1}+x_{2}$, and in general, we are looking at $\mathbb{R}^{n}\mapsto\mathbb{R}$ specified by $(x_{1},x_{2},\ldots,x_{n})\mapsto\sum_{i=1}^{n}x_{i}$. Squaring both sides, it is clear that

On the other hand, since

the addition of edge lengths can result in an expansion bounded by the size of the sum.

For a rooted tree $T$, let $0pt(T)$ denote the length of the longest path from a leaf to the root.

Using Proposition 4.8, it is straighforward to deduce the next two theorems that provide the theoretical support for the use of tree dimensionality reduction. The following theorem is an immediate consequence of Proposition 4.8, choosing the path realizing the distance between $T$ and $T^{\prime}$.

Let $A$ and $B$ be subsets of $\mathcal{E}_{m}$ such that each item of $A$ and $B$ has a label in $L\subset\mathcal{P}(\{1,\ldots,m\})$ (where $\mathcal{P}(-)$ denotes the power set of $\{1,\ldots,m\}$). Then we can define a matching distance as

Without assuming such a labelling, we define the matching distance between $A$ and $B$ to be

where $\phi$ varies over all bijections $A\to B$.

The following is now also immediate from Proposition 4.8.

We conclude the discussion by describing some computational results on simulated data that illustrates the divergence between $d_{M}$ and $d_{\mathcal{E}_{m}}(T,T^{\prime})$; see Figure 7. To demonstrate the stability of the tree projection operation, and to empirically test its approach of the distance bound, we constructed a panel of pairs of $m$-dimensional trees, $m>10$. The distances between the pairs of trees were computed and compared against the distributions of distances induced by the $k$-dimensional projection operator, $\Psi_{k}$ with $k\in{3,4,5,6}$. The distances between elements of the projections can exceed the original inter-phylogeny distance ($\epsilon$), but rarely approach the upper bound. To further characterize the behavior of the distribution of subtree distances as a function of $k$, we compared the distributions of projected distances ranging from $k=3$ to $k=m-1$ for a fixed pair of $m$-dimensional phylogenies. We see that the median approaches the true $m$-dimensional $d_{BHV}(T,T^{\prime})$ with larger values of $k$, while the variance appears to decrease monotonically.

Broadly speaking, the various tree dimensionality reduction operators transform questions about comparison of trees or analysis of finite sets of trees to questions about comparisons and analysis of sets of clouds of trees. This has several advantages. First, when working with $\mathbb{P}\Sigma_{m}$, the resulting clouds live in $\Sigma_{k}$, and as discussed above analysis in the non-projectivized space can be simpler. Second, when projecting to $\Sigma_{3}$ or $\Sigma_{4}$, both visualization and analysis is easier (particularly in $\Sigma_{3}$, since that has a Euclidean metric). For example, in order to perform supervised classification on a labelled set of trees $X$, we can simultaneously solve classification problems in $\Psi^{\prime}_{k}(X)$ and use majority voting in order to assign labels to new trees. Another possibility is to perform hierarchical clustering on the clouds in $\Psi^{\prime}_{k}(X)$, resulting in another tree, and use the distance in tree space as a test statistic to discriminate between clouds. In Section 6 below we describe an application that involves producing a predictor for influenza vaccine effectiveness using the variance of the distribution produced by $\Psi$.

Chronic lymphocytic leukemia (CLL) is the most common leukemia in adults, primarily affecting the elderly population (median age at diagnosis is 70) [51]. CLL is a proliferative disorder of B-lymphoctyes characterized by a steady accumulation of clonal, non-functional B-cells. Treatment strategies vary greatly given the heterogeneity in disease course, ranging from watchful waiting, to localized radiation, to systemic chemotherapy. The fact that CLL is a relatively indolent malignancy makes it an excellent model for studying clonal evolution under different therapeutic strategies [56].

A recent genomic study [31] performed whole exome sequencing on 160 CLL cases covering the spectrum of clinical courses the disease can take. This data established a space of recurrent alterations, which was then used to genotype 18 patients for whom two time points were available. Of these 18 patients, 10 of 12 treated with chemotherapy underwent clonal evolution compared to only 1 of 6 receiving no treatment according to the authors of the study. We combine the 18 patients of [31] with those of a similar study [47] wherein 3 CLL patients received chemotherapy and were sequenced at multiple time points. The multiple time points for the 3 patients studied in [47] are decomposed into all combinatorial triplets. Therefore, there are a total of 18 phylogenetic trees inferred from [31] and 12 phylogenetic trees inferred from [47]. In Figure 8, we map this data to $\mathbb{P}\Sigma_{3}$ and color based on treatment status.

From a clinical standpoint, a central question is whether treatment promotes evolution of the cancer and whether there is strong evidence for avoiding cytotoxic therapy in patient management. Quite clearly the distribution of 6 untreated patients resembles a pattern (in green) forms a tight cluster, indicating that tumors from patients who did not received therapy are stable genetically, sharing most of the mutations. However, in red are represented the histories of tumors of 15 patients under therapy, presenting some mutations at different times that are not shared across different samples. The ratio between the number of mutations that are exclusive in the early branch versus the ones that are share with other phases is represented in right hand side of Figure 8. A number close to zero indicates a genetically stable tumor.

To assess how different are the clonal histories of tumors from untreated vs treated patients, we studied the distance between the centroids of the two populations regarded as points in $\mathbb{P}\Sigma_{3}$. The 95% CI for the distance between the centroids of the treated / untreated groups is (0.15, 0.36), under 1000-fold bootstrap resampling. The analogous intervals for untreated / untreated and treated / treated are (0.01, 0.14) and (0.02, 0.16) respectively. This analysis shows that the centroids of these clinically distinct sets of patients are well-resolved, supporting the idea that untreated tumors are more stable than treated ones, where district mutations can appear along the evolution of the tumor.

As another application, we examined low grade gliomas (LGG), a set of tumors of the central nervous system most often involving astrocytes or oligodendrocytes. They are distinguised from high grade gliomas (III, IV), such as glioblastoma multiforme, by the absence of anaplasia and have a more favorable prognosis. Surgery alone is not considered curative for LGG and patients are typically treated with adjuvant radiation therapy, chemotherapy, or both. If a patient relapses, the tumor may be observed to have a higher grade at that time. Johnson et al. studied a cohort of 23 LGG patients who relapsed, many of whom were treated with the chemotherapeutic agent temozolamide (TMZ) [29]. Whole exome sequencing was performed on tumor tissue at diagnosis and at relapse in an effort to characterize the evolution of recurrent glioma.

TMZ is an alkylating agent that directly damages the cellular genome, and accordingly we see that the subset of patients that were treated with TMZ show a greater acquisition of relapse-specific mutations. After projecting the trees to $\mathbb{P}\Sigma_{3}$, we find the 95% confidence interval for distance between centroids of the treated / untreated groups is (0.31, 0.48), under 1000-fold bootstrap resampling. The analogous intervals for untreated / untreated and treated / treated are (0.02, 0.13) and (0.02, 0.16) respectively. TMZ treatment status defines two statistically well-resolved sub-populations of patients with respect to the shape of their evolutionary behavior, an observation recently reinforced by the larger study of [55]. Also of interest is the apparent correlation between the geometry of the phylogenetic tree and the histologic grade at relapse. A 1-nearest neighbor classifier of this trinary observable (grade II, III, or IV at relapse) yields 85% accuracy under two-fold cross-validation, and the accuracy does not improve with larger values of $k$. The tesselation of the space associated to a 1-nearest neighbor classifier is known as a Voronoi diagram, and we have colored the cells in accordance with the grade at relapse.

This example shows that in the simple case of trees with three leaves, evolutionary tumor histories are visibly different under different therapeutic regimes. Moreover, we see that evolutionary trajectory can be associated with prognosis, as measured by grade at relapse.

The recent development of single cell transcriptomics and genomics is providing an opportunity to study the role of clonal heterogeneity in tumors [34, 16, 40] and to identify small, previously uncharacterized cell populations [24]. The single cell approach to studying complex populations brings with it new challenges associated with the large number of sampled genomes. Another rapidly maturing technology in the modeling of tumor evolution is that of patient-derived xenografts. Patient-derived xenografts (PDX) are generated by transplanting tumor tissue into immunodeficient mice, serially engrafting in new mice as each host animal expires. This provides an in vivo platform for drug screening as well as longitudinal monitoring of tumor adaptation and clonal dynamics.

We take advantage of recently published PDX data from breast cancer patients, where single-nucleus deep-sequencing was performed across a lineage of host animals engrafted with a primary lesion of triple-negative breast cancer [16]. The data are comprised of normal tissue, a sample from the primary tumor, and three subsequent mouse passages. Somatic mutation calls revealed 55 informative sites of substitution in this lineage, and these variants were assigned to eight distinct cellular populations based on bulk sequencing.

Only bulk sequencing data were available from the primary tumor and matched normal tissue, while the three xenograft passages were sequenced at single-nucleus resolution. For each cellular fraction we randomly sampled a single nucleus from the first, second, and fourth PDX passages (27, 36, 27 nuclei respectively were available). This preparation of the data implies five sequential time points along the tumor’s history: benign germline genome, genotype at diagnosis of primary tumor, genotype at first PDX, genotype at second PDX, and genotype at fourth PDX. The combinatorial possibilities in the latter three time points give rise to a large forest of phylogenetic trees.

We examined the difference in heterogeneity between phylogenetic trees constructed on the first four time points and those constructed on last four time points. The distributions trees from both time windows are visualized in Figure 10 as point clouds in $\mathbb{P}\Sigma_{4}$. Consistent linear evolution is seen from primary tumor through the first two xenograft passages, however we observe significant heterogeneity of tumor clones upon the fourth mouse passage. The first time window (purple) is completely contained within the topology corresponding to linear evolution, unlike the second (gold) which is centered on the origin and extends into all three possible topologies. The point cloud for the second time window displays a higher standard deviation than the first (10.49 vs. 8.69), and its centroid is essentially a star tree. Centroids and variances are computed in $\Sigma_{4}$, prior to rescaling of branch lengths. The high degree of genotypic heterogeneity giving rise to the second time window distribution is suggestive of a clonal replacement event between the time points of Xenograft 2 (X2) and Xenograft 4 (X4). Many of the prevalent alterations before X4 disappear during the final passage, and many new mutations rise to dominance. These results raise interesting questions about the long-term fidelity of PDX vehicles to the genetics of their ancestral primary tumors, which theoretically they serve to mimic.

The authors gratefully acknowledge the constructive feedback of Gillian Grindstaff, Melissa McGuirl, and Daniel Rosenbloom. This work was supported by a TL1 personalized medicine fellowship (5TL1TR000082) and NIH grants (R01 CA179044, R01CA185486, R01GM117591, U54 CA193313).