Background#
This page defines the problem ForestPD solves and outlines the algorithm. The manuscript (see the citation on the front page) has the definitions, proofs and experiments.
Phylogenetic diversity#
One tree#
Let \(X\) be a finite set of taxa. A phylogenetic tree \(\mathcal{T} = (T, \omega)\) on \(X\) consists of a tree \(T\) whose leaves are labelled bijectively with the taxa in \(X\), and an edge weight function \(\omega \colon E(T) \to \mathbb{Q}_{\geq 0}\) (the branch lengths). An internal vertex of degree larger than 3 is a polytomy, and a tree without polytomies is binary. For a set \(A \subseteq X\), \(E_{\mathcal{T}}(A)\) is the edge set of the spanning tree of \(A\), the smallest subtree of \(\mathcal{T}\) that contains the leaves in \(A\). The phylogenetic diversity (PD) of \(A\) is [Faith1992]
the total weight of the edges on paths between taxa in \(A\). Closely related taxa share most of their branches, so protecting both adds little over protecting one. PD therefore favours sets of complementary taxa, which is why conservation planning maximizes the PD of a set rather than ranking taxa individually.
By default PD is unrooted, as defined above, so a single taxon has PD 0. With rooted=True
it is rooted PD: the branches from the chosen taxa to the root, plus the root’s own branch,
count (see Rooted PD).
On one tree, a set of \(m\) taxa with the largest PD can be found greedily [Steel2005] [Pardi2005], and tools such as PDA solve this for large trees within seconds [Minh2006].
Many trees#
A phylogenetic analysis usually gives several trees: a Bayesian posterior sample, bootstrap replicates, or gene trees that differ from each other and from the species tree, for example through incomplete lineage sorting. ForestPD maximizes the weighted sum of PD over a collection \(\mathcal{T} = \{\mathcal{T}_1, \dots, \mathcal{T}_t\}\) of trees on \(X\), instead of the PD on a single or summary tree:
The tree weights \(w(\mathcal{T}_i) > 0\) can be posterior probabilities (the objective is then the expected PD under phylogenetic uncertainty), \(1/t\) (the average PD over the trees), or any other weighting, such as alignment lengths for gene trees.
For example, the two trees below have the same four taxa but different topologies. On \(\mathcal{T}_1\) every pair across the middle edge, such as \(\{A, C\}\), has the largest PD, 4. In \(\mathcal{T}_2\) the taxa \(A\) and \(C\) are siblings, and the same pair has PD 2. Only \(\{A, D\}\) and \(\{B, C\}\) have PD 4 in both trees, for a total of 8:
Two trees on taxa A–D, with branch lengths. The highlighted edges connect \(A\) and \(D\): PD 4 in each tree, 8 in total. \(\{A, C\}\) has PD 4 in \(\mathcal{T}_1\) but only 2 in \(\mathcal{T}_2\).#
>>> import forestpd as fpd
>>> trees = ["(A:1,B:1,(C:1,D:1):2);", "(A:1,C:1,(B:1,D:1):2);"]
>>> instance = fpd.MaxMPDInstance(trees)
>>> instance.evaluate_pd({"A", "C"}), instance.evaluate_pd({"A", "D"})
(6.0, 8.0)
>>> instance.optimize(budget=2).value
8.0
Costs and a budget#
Each taxon \(x\) has an integer cost \(c(x) \geq 1\), and the protected set \(S\) must fit a budget \(B\): \(c(S) \leq B\), where \(c(S) = \sum_{x \in S} c(x)\). With unit costs (the default) the budget is the number of taxa.
Extinction risk#
Each taxon \(x\) has a survival probability \(a(x)\): the probability that it survives without protection. For a protected set \(S\), let \(a_S(x) = 1\) if \(x \in S\) and \(a_S(x) = a(x)\) otherwise. Taxa survive independently, so the probability that exactly the taxa in \(A\) survive is
and the expected weighted PD of the survivors is [Weitzman1998]
An edge contributes if a taxon survives on each of its sides. ForestPD therefore computes the expectation per edge: an edge \(e\) that splits \(X\) into \(X_e\) and \(\overline{X}_e\) contributes \(\omega(e)\) times \(\big(1 - \prod_{x \in X_e} (1 - a_S(x))\big) \big(1 - \prod_{x \in \overline{X}_e} (1 - a_S(x))\big)\). With every \(a(x) = 0\) (the default), the objective is the weighted PD of \(S\).
MaxMPD#
Maximize Multi-Tree Phylogenetic Diversity (MaxMPD)
Input: a set of taxa \(X\); a collection \(\mathcal{T} = \{\mathcal{T}_1, \dots, \mathcal{T}_t\}\) of phylogenetic trees on \(X\); tree weights \(w \colon \mathcal{T} \to \mathbb{Q}_{> 0}\); taxon costs \(c \colon X \to \mathbb{N}_{> 0}\); survival probabilities \(a \colon X \to [0, 1) \cap \mathbb{Q}\); and a budget \(B \in \mathbb{N}_{> 0}\).
Output: a set \(S \subseteq X\) with \(c(S) \leq B\) that maximizes \(\sum_{A \subseteq X} a_S(A) \sum_{i=1}^{t} w(\mathcal{T}_i) \cdot \mathrm{PD}_{\mathcal{T}_i}(A)\).
We write \(n = |X|\) for the number of taxa and \(t = |\mathcal{T}|\) for the number of
trees. The trees may be rooted or unrooted (with rooted=True) and binary or not, but every tree
must contain every taxon. MaxMPD is NP-hard already for three trees. ForestPD’s algorithm is fast
when the trees differ only locally, as trees inferred from the same data typically do.
Algorithm#
ForestPD solves MaxMPD exactly, in three stages.
Cluster reduction#
Wherever every tree has the same split of the taxa (an edge separating the same two groups), the instance is cut at that split into smaller instances, one per cluster, with a placeholder leaf for the rest of the trees. The clusters form a tree and are solved from the leaves up: each cluster passes its optimal value per budget to its parent. Trees that share most of their structure give many small clusters.
Agreement forests#
Within a cluster, the trees still differ. An agreement forest partitions the cluster’s leaves
into \(k\) components such that every component has the same topology in every tree, and in
each tree the components occupy disjoint parts. Its size \(k\) measures how much the trees
differ in that cluster: \(k = 1\) means they agree completely. Finding a smallest agreement
forest is NP-hard, but every agreement forest gives the exact optimum; a smaller one makes the next
stage faster. ForestPD’s default is a fast heuristic (HeuristicMAF); with
the optional maffe package it can use minimum agreement
forests (see Agreement-forest stages).
Dynamic programming over signatures#
For each cluster and each signature, the subset of components that contain a protected taxon, the objective splits into a constant plus independent contributions of the components. Each component’s contribution is computed for every budget by a dynamic program along its trees, and the budgets are divided over the components by max-plus convolution.
Running time#
With \(2^k\) signatures per cluster, the running time is
for budget \(B\) and total cost \(C = c(X)\), where \(k\) is the largest agreement-forest size over the clusters. The \(\min(B, C - B)\) comes from two kinds of tables: ones indexed by the budget spent (up to \(B\)), and complementary ones indexed by the budget left unspent (up to \(C - B\)). ForestPD computes whichever is smaller for the requested budgets, or both. In practice:
Budget curves. The tables cover every budget up to \(B\), so one run gives the optimum for every smaller budget as well (see Several budgets).
Hard instances have a large \(k\) in some cluster, which occurs where many taxa are placed differently by different trees, and, to a lesser extent, a large budget in cost units (fine-grained costs). The numbers of trees and taxa enter only linearly.
structure()reports \(k\) before solving.
References#
D. P. Faith. Conservation evaluation and phylogenetic diversity. Biological Conservation 61(1):1–10, 1992.
M. Steel. Phylogenetic diversity and the greedy algorithm. Systematic Biology 54(4):527–529, 2005.
F. Pardi and N. Goldman. Species choice for comparative genomics: being greedy works. PLoS Genetics 1(6):e71, 2005.
B. Q. Minh, S. Klaere, and A. von Haeseler. Phylogenetic diversity within seconds. Systematic Biology 55(5):769–773, 2006.
Weitzman. The Noah’s Ark problem. Econometrica 66(6):1279–1298, 1998.