Can you trust your reconstructed lineage tree? A homoplasy-based approach for irreversible evolution

Can you trust your reconstructed lineage tree? A homoplasy-based approach for irreversible evolution | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Can you trust your reconstructed lineage tree? A homoplasy-based approach for irreversible evolution View ORCID Profile Pini Zilber , View ORCID Profile Sebastian Prillo , View ORCID Profile Nir Yosef , View ORCID Profile Boaz Nadler doi: https://doi.org/10.1101/2025.07.27.667007 Pini Zilber * Department of Computer Science and Applied Mathematics, Weizmann Institute of Science , Israel Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Pini Zilber Sebastian Prillo † Department of Electrical Engineering and Computer Sciences, University of California , Berkeley, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sebastian Prillo Nir Yosef † Department of Electrical Engineering and Computer Sciences, University of California , Berkeley, USA ‡ Department of Systems Immunology, Weizmann Institute of Science , Israel Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Nir Yosef Boaz Nadler * Department of Computer Science and Applied Mathematics, Weizmann Institute of Science , Israel Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Boaz Nadler Abstract Full Text Info/History Metrics Preview PDF Abstract Phylogeny inference is a fundamental problem in computational biology, with many proposed algorithms. Emerging techniques that couple single-cell genomics with Cas9-based genome editing open the way for indepth analysis of cell phylogenies that underlie processes of clonal expansion, selection and diversification, from embryogenesis to cancer. A key distinguishing feature of cell lineage analysis with these techniques is the non-modifiability of Cas9-induced mutations, which motivates revisiting questions in phylogenetics. In this work, we ask one such fundamental question: is it possible to assess the reliability of an inferred lineage tree, even though we do not know its underlying ground truth? We present a homoplasy-based approach for this question that leverages the non-modifiability property. We show via simulations that under a broad range of settings, our method can effectively distinguish accurate reconstructions out of a pool of candidate solutions. Importantly, our homoplasy-based score is substantially more powerful than the commonly used parsimony score - a result that we back by both empirical and theoretical analysis. The computation of the homoplasy score is simple and scalable, thus opening the way for more rigorous analysis of cell lineages. 1 Introduction Reconstruction of phylogenies has long stood as a cornerstone of evolutionary biology, with numerous methods developed over the years [ Hen99 , NK00 , Hal04 , WL11 ]. In the classical problem, the objective is to reconstruct the tree of life. The leaves of the tree represent extant species, and the internal nodes provide a way to reason about their latent ancestors far back in evolutionary history. Phylogenetics has also been used to formalize the analysis of sequences of nucleic or amino acids and shed light on the evolutionary forces that shaped them. The process of natural evolution of sequences, captured by classical models such as Jukes-Cantor, reflects the chances of mutations to occur and is assumed to be modifiable . Namely, the occurrence of a mutation in a site does not necessarily nullify the chances of that site being mutated again. The absence of this property is one of the major distinguishing features of an emerging branch of phylogenetics, where trees depict the evolution of cellular clones, with leaves corresponding to cells and internal nodes reflecting their sub-clonal relationships [ GGH+21 ]. While cell lineages can be inferred based on somatic mutations or other natural and modifiable forms of heritable variation, the most powerful techniques (in terms of numbers of cells and the resolution of lineage trees) are based on inducible mutations at synthetic DNA elements [ KVO18 , MG19 , YJN+22 ]. In these applications, the heritable information that guides the retrospective inference of lineages is accrued through genome editing (e.g. with CRISPR/Cas9) of synthetic “recorder sites” that are integrated into the genome (usually at the 3’ end of synthetic genes). Since CRISPR/Cas9 often has a much lower affinity to the edited sites, it is unlikely for additional mutations to occur in already-mutated sites, making this process non-modifiable [ SSCR23 ]. Several studies investigated this new setting and proposed suitable algorithms to reconstruct cell lineages based on mutation profiles that accumulated in each cell [ JKQ + 20 , FDM + 21 , SS22 , WZKY23 , SSCR23 ] (see also [ GGH + 21 ] for a broad performance comparison). One strategy relied on the long-studied Camin-Sokal model [ CS65 ]. This model assumes the less restrictive scenario of irreversibility , which permits mutated sites to undergo further mutations, but not to revert to an ancestral state. This is a notable difference from non-modifiability, except for the special case of a two-state evolution, where the two become identical. We note that in practical applications, other algorithms, including those designed for the classical problem such as neighbor-joining [ SN87 ], are often effective [ SL19 , JKQ + 20 , PAW + 25]. This multitude of strategies and algorithms for the inference problem raises the following question: given several reconstructions of the cell lineage tree that are all based on the same (non-modifiable) data, which reconstruction is more accurate? (i.e. more similar to the true one). More fundamentally, can we assess whether any individual reconstructed topology is accurate or not? This latter objective extends beyond simply ranking candidate topologies since it also aims to determine if any of the candidate topologies, e.g. the top-ranked one, are close enough to the true topology and can be used for downstream analysis. A common strategy to evaluate candidate solutions assumes that the true topology is known and compares directly to it (for example, using the Robinson-Foulds distance or the frequency of correct triplets [ JKQ + 20 , GGH + 21 ]). While this strategy is naturally restricted, other approaches have been proposed that do not assume any knowledge of the ground truth, including comparison of parsimony (the minimal number of ancestral mutation events that explain the data), or likelihood (given an evolutionary model of mutation accrual). These approaches, however, are often used only for ranking solutions and not for making “absolute” decisions of whether or not a given tree is similar to the latent ground truth. Furthermore, their ability to effectively distinguish accurate from inaccurate solutions is unclear. In this work, we present the pairwise homoplasy score (PHS) - a simple and theoretically grounded approach to address these problems. Leveraging the non-modifiable nature of the evolution process, our algorithm calculates a sequence of tail probabilities - one for each pair of leaves in the reconstructed tree - and combines them into a single score. We demonstrate that compared to parsimony and likelihood, the PHS approach is substantially more powerful in distinguishing accurate tree reconstructions out of a pool of candidate solutions. Furthermore, we derive a theoretical analysis that provides insight into the advantage of PHS over parsimony. The procedure for calculating the PHS has polynomial-time complexity, and can be easily run on trees with thousands of leaves. 1 We expect it to become an important part of future pipelines for phylogenetic analysis in non-modifiable systems, and particularly CRISPR/Cas9-based tracing of cell lineages. 2 Problem Setup 2.1 CRISPR-Cas9 non-modifiable mutation model We consider the following model for a CRISPR-Cas9 cell lineage tracing experiment, similar to [ WZKY23 ]. A typical experiment starts with a single progenitor cell that has k target sites along its genome, which are all unedited (i.e., yet unmutated). As the experiment ensues, the resulting progeny grows through cell divisions, with possible divergence from neutral growth due to advantageous properties of some clades. This evolutionary process can be described by a binary tree, where nodes correspond to cells. Each node u has a birth time τ u and a sequence s ( u ) of k characters, which are the Cas9-induced mutations at its k sites. The root node of the tree, denoted by r , corresponds to the progenitor cell. It has a birth time τ r = 0 and an unmutated sequence s r = (0, 0, …, 0). The path length τ u,v between a node u and a descendant node v is defined as the difference between their birth times, τ u,v = τ v − τ u . Following [ WZKY23 ], we assume that characters evolve independently along the tree. Specifically, each unmutated character i ∈ [ k ] (namely with s i = 0) evolves according to a Poisson process with a mutation rate λ . For simplicity, we assume this rate λ is fixed for all character locations in the tree; some extensions are discussed in Section 4 . Once a mutation occurs, the mutated state is drawn from a probability distribution over a set of m possible states 1, …, m where m ≥ 2. The evolution is non-modifiable : once a character is mutated, it cannot mutate again and its state remains fixed. Furthermore, all mutations of a cell are passed on to all its descendants. At the end of the experiment, the sequences at a subset of n cells are observed. Without loss of generality, we set the time of the end of the experiment at τ = 1. We denote the n×k matrix of the sequences at the n cells by 𝒮 n . These n cells induce an underlying ground-truth binary tree, denoted 𝒯 GT , which describes their clonal history. For future use, we denote the set of all node sequences in this tree by 𝒮 GT ={ s ( u ) : u is a node in 𝒯 GT }. By definition, 𝒮 n is the subset of 𝒮 GT that corresponds to its terminal leaves. For simplicity, we assume that 𝒮 n is fully observed. Settings with missing data, due to either stochastic or heritable missingness, are addressed in Appendix B. There are two key quantities that characterize the non-modifiable process described above. One is the probability ρ that an observed character at the end of the experiment is mutated. It is given by The second quantity is the collision probability, denoted by q . This is the probability that two independent mutation events result in the same mutated state, and it is given by Since there are m ≥ 2 possible mutated states, the collision probability satisfies q < 1. All the parameters characterizing the tree topology and the generative process for the sequences are summarized in Table 1 . We conclude this subsection with the following definition that will be used extensively in the manuscript. View this table: View inline View popup Download powerpoint Table 1: Model parameters. Definition 1 (full tree). A full tree 𝒯 𝒮 is a structure that consists of both a tree topology 𝒯, branch lengths (elapsed time between divisions), and a collection of the character sequences 𝒮 at all its tree nodes. 2.2 Distinguishing between accurate and inaccurate reconstructed trees Given the n observed sequences 𝒮 n , key problems are to reconstruct either the underlying tree topology 𝒯 GT or the full tree 𝒯 GT 𝒮 GT . Several reconstruction algorithms were developed for non-modifiable CRISPR-Cas9 evolution models, for example [ JKQ + 20 , FDM + 21 , SS22 , WZKY23 , SSCR23 ]. The celebrated Neighbor-Joining method [ SN87 ], originally developed for reversible evolution models, was shown to perform well also in non-modifiable settings [ SL19 , JKQ + 20 , PAW + 25 ]. Given the input sequences 𝒮 n , different reconstruction algorithms often output different tree topologies. Some of these reconstructed trees may be accurate, whereas others may be quite distant from the ground-truth tree. This may be due to the limited sequence length k and the fact that various methods take different optimization approaches and approximations to the NP-hard problem of tree reconstruction. The problem at the focus of our work is to develop a method to distinguish whether a given reconstructed tree is accurate, or inaccurate and far from the ground-truth tree. Crucially, this task needs to be accomplished with the ground-truth tree remaining unknown. In the phylogenetics literature, several works considered a different, though related problem of proposing meaningful distances between a reconstructed tree 𝒯 and a known ground truth 𝒯 GT , such as the Robinson-Foulds distance [ RF81 ] and the triplets score [ JKQ + 20 ]. These methods are useful in simulation studies, but not for assessing tree reconstructions in practical settings, where 𝒯 GT is unknown. Two common measures that do not require knowledge of 𝒯 GT are the parsimony and the likelihood of a tree. In the parsimony approach, a simple polynomial-time procedure is applied to a reconstructed topology 𝒯 to estimate the minimal number of ancestral mutations that give rise to the observed leaf states [ Gus97 ]; see Appendix C. In contrast, the likelihood approach considers all possible combinations of ancestral states. Under a Markov assumption, this score can be calculated efficiently, while further requiring a model of mutation accrual and knowledge of edge lengths [ Gus97 ]. According to the maximum parsimony principle, trees with fewer number of mutations M (𝒯 𝒮) are considered closer to the ground-truth one [ Fit71 , Alb05 ]. Similarly, trees with higher likelihood values are considered more accurate. Likelihood is more statistically grounded than parsimony: Under the assumption of a generative model of evolution, trees with maximum likelihood L (𝒯 𝒮) are asymptotically consistent [ Fel81 , War18 ]. In principle, a set of reconstructed trees may be ordered by their parsimony or likelihood values. However, it is not possible to tell from the parsimony and likelihood values which trees are accurate and which are not. In particular, if all trees are far from the ground truth, this ordering may be meaningless. We empirically illustrate these issues in sections 3 and 5 . A fundamental challenge is thus to devise a score that is able to detect if a given tree is accurate or not, for example indicate if its normalized RF distance from the unknown ground-truth tree is smaller than some prescribed value ϵ . We mathematically formulate this challenge as follows: Given a tree 𝒯 𝒮, a distance measure d and a threshold ϵ ∈ (0, 1), we consider the following hypothesis testing problem, We emphasize that the main obstacle is that 𝒯 GT 𝒮 GT , and thus d (𝒯 𝒮, 𝒯 GT 𝒮 GT ), are unknown. Note that in (3), the null hypothesis is the desirable one, in which the candidate tree is close to the ground-truth one. In this paper, we consider the following four distance functions: the normalized RF distance the triplets distance a parsimony-based distance and (iv) a likelihood-based distance, The max operator in Eqs. (6) and (7) ensures that these distances are non-negative. This is required as, given observed data 𝒮 n , the ground-truth tree 𝒯 GT 𝒮 GT might not be the one with maximum parsimony or highest likelihood. Our main contribution is an approach to resolve (3), which is applicable and powerful for a wide choice of distance functions between trees and of cutoff values ϵ . Specifically, we propose a homoplasy-based test statistic that quantifies the consistency of the candidate tree with respect to the observed sequences 𝒮 n and the non-modifiability of the CRISPR-Cas9 model. As we show, our approach can distinguish accurate from inaccurate trees, under all four distance functions d RF , d tri , d P and d L . In Section 5 we present simulation results for d RF and d P ; simulations for d tri and d L appear in the appendix. In addition, in Section 6 we provide theoretical support for our approach. Remark 1. There is a fundamental difference between the problem defined in Eq. (3) and classical statistical hypothesis testing. The latter is often formulated as a decision problem, whether observed data was generated from an assumed statistical model (the null) or an alternative one. In this work, in contrast, we assume that the CRISPR-Cas9 statistical model is correct. Instead, we test the consistency of the “data” 𝒯 𝒮 with the CRISPR-Cas9 model. Here, the “data” includes both the observed sequences and the reconstructed tree and its inner sequences. 3 Parsimony and Likelihood Baselines To motivate our approach, we first demonstrate via simulations that given a reconstructed tree 𝒯 𝒮 corresponding to an unknown ground truth 𝒯 GT 𝒮 GT , existing measures are unable to accurately distinguish between the null and the alternative hypotheses in (3). In what follows, we consider two measures based on parsimony and one based on the likelihood of the reconstructed tree. A first approach is to use the parsimony score of the reconstructed tree, M (𝒯 𝒮). Specifically, for a suitably chosen cutoff value t > 0, the reconstructed tree 𝒯 𝒮 is rejected if M (𝒯 𝒮) > t and accepted otherwise. In general, this threshold t may depend on the model parameters, as well as on the specific distance function in (3), which may be different from the parsimony distance d P . A second measure is based on the parsimony probability distribution conditional on the specific reconstructed tree topology 𝒯, denoted as ℙ[ M | 𝒯]. As described in Appendix D, this distribution can be computed analytically. Given the reconstructed tree 𝒯 and its reconstructed sequences 𝒮, we calculate the tail probability p (𝒯 𝒮) = ∑ M≥M ( 𝒯 𝒮 ) ℙ[ M | 𝒯]. A reconstructed tree 𝒯 𝒮 is rejected if p (𝒯 𝒮) < t for some suitable threshold 0 < t < 1. Finally, a third measure is based on the likelihood L (𝒯 𝒮). Specifically, it rejects the reconstructed tree 𝒯 𝒮 if its log-likelihood, log L (𝒯 𝒮), is below a cutoff value t . Next, we illustrate via several simulations the limitations of these three approaches in resolving the hypotheses in Eq. (3) . To this end, we generated 1000 ground-truth trees and their inner sequences, 𝒯 GT 𝒮 GT . Similar to other studies [ JKQ + 20 , FDM + 21 , SS22 , SSCR23 ], each ground-truth tree topology was randomly generated according to a birth-death process. The rates of birth and death are not homogeneous throughout the tree but instead depend on a heritable, sub-clonal fitness level; see Appendix J. Given a tree topology, its inner sequences were generated by the non-modifiable CRISPR-Cas9 mutation process described in Section 2.1 . For each 𝒯 GT 𝒮 GT , we run several reconstruction algorithms, listed in Section 5 , with the observed data S n given as their input. In our first simulation, we considered a tree as inaccurate if M (𝒯 𝒮) > 1.01 · M (𝒯 GT 𝒮 GT ). This corresponds to the hypothesis testing problem in Eq. (3) with distance function d = d P and ϵ = 0.01; results for other distances d and thresholds ϵ appear in Section 5 and Appendix K. As shown in panels (a-c) of Fig. 1 , for each of these three baseline measures there is substantial overlap between their values on accurate and inaccurate reconstructed trees. A similar overlap exists even between ground-truth and inaccurate reconstructed trees. This substantial overlap implies that regardless of the chosen threshold t , these measures cannot accurately separate ground-truth trees, or accurate trees, from inaccurate ones. Note that even when accuracy is defined by the parsimony -based distance d P , both the parsimony score and its tail-probability measure are ineffective in distinguishing between accurate and inaccurate reconstructed trees. Download figure Open in new tab Figure 1: Kernel density estimates for four accuracy measures: parsimony (upper left); parsimony tail-probability (upper right); likelihood (bottom left); and our proposed cPHS (bottom right). The parsimony tail probability and cPHS test statistics are log-transformed for clarity of view. The hypotheses ℋ 0 and ℋ 1 are defined by (3) with d = d P (6) and ϵ = 0.01. The simulations were performed using Cassiopeia [ JKQ + 20 ] as follows: binary tres with 10 6 terminal cells were generated, from which n = 10 3 observed cells were subsampled (subsampling ratio of 10 −3 ). The sequence length is k = 50, each site has m = 50 possible mutations with probabilities q i drawn from an exponential distribution, and λ = −log 2 so that there is ρ = 50% probability to observe a mutation at a leaf. These parameter values are similar to those of experimental settings, see [ JKQ + 20 ]. In contrast, Fig. 1(d) shows that our proposed PHS-based test statistic, described in the next section, achieves a high separation between candidate trees with low (good) and high (bad) parsimony. Furthermore, it achieves a perfect separation between ground-truth trees and candidate trees with high parsimony. In Section 5 and Appendix K we show that a similar conclusion holds for the other distances d RF , d tri , and d L . Specifically, for each of these other distance functions, there is significant overlap in the distribution based on the parsimony or likelihood measures. In contrast, our PHS-based test statistic achieves a high separability. In Section 6 we provide theoretical support for the empirical advantage of our homoplasy-based approach over parsimony. Specifically, we prove that to detect certain incorrect topologies, parsimony requires exponentially longer sequences (in tree depth) compared to our PHS approach. 4 A Pairwise Homoplasy Approach To describe our approach for the hypothesis testing in Eq. (3) , we first recall the definition of two classical concepts in tree phylogeny: latest (i.e. lowest or most recent) common ancestor (LCA) and homoplasy. Definition 2 (LCA). The latest common ancestor of a pair of leaves u, v with respect to a tree topology 𝒯 is the most recent node (with latest birth time τ ) whose descendants include both u and v , denoted LCA( u, v ). Definition 3 (homoplasy and PHS). For a pair of leaves u and v , there is a homoplasy in their i -th character with respect to a full tree 𝒯 𝒮 if both leaves have the same mutation, , but their latest common ancestor was unmutated, . We denote the indicator of such a homoplasy in the i -th character of u and v by and by phs( u, v ) the pairwise homoplasy score for the pair of leaves u, v with respect to a full tree 𝒯 𝒮, defined as their mean PHS over the k characters, For example, consider a full tree 𝒯 𝒮 with s ( u ) = 10221, s ( v ) = 13321, and s (LCA(u,v)) = 00020. Then phs 1 ( u, v ) = phs 5 ( u, v ) = 1, while phs 2 ( u, v ) = phs 3 ( u, v ) = phs 4 ( u, v ) = 0, and thus phs( u, v ) = 2 / 5. Note that to compute the PHS, the unobserved sequence s ( w ) at the inner node w = LCA( u, v ) needs tobe known. In cases where a reconstruction algorithm provides us only with a tree topology 𝒯 but not with its estimated inner sequences 𝒮, we first impute the ancestral states using a simple procedure of post-order traversal (Appendix A.2). Given two sister nodes, the state of their parent can in most cases be inferred without ambiguity due to non-modifiability. The one exception is when the two sister nodes are mutated and have the same state. In this case, we make the parsimonious choice and set the parent cell to the same state as its child nodes. 4.1 PHS probability distribution Let 𝒯 GT be a fixed ground-truth topology, and suppose its sequences 𝒮 GT are generated according to the non-modifiable mutation model of Section 2.1 . As the sequences are random, so are the PHS values of all pairs of tree leaves. Our test statistic is based on the tail probabilities of these values. To construct our test statistic, we shall make use of the following lemma, which describes the PHS probability distribution. Its proof is in Appendix E. Lemma 1. Let u and v be a pair of leaves in 𝒯 GT . Denote their LCA by w = LCA ( u, v ) and its birth time by τ w ∈ [0, 1). Let and , where λ is the mutation rate. Then, for sequences of length k generated as in Section 2.1 with collision probability q (see Eq. (2) ) , Lemma 1 reveals two appealing properties of the PHS distribution. First, in terms of the model parameters, it depends on the mutation rate λ , sequence length k and the collision probability q , but not on the individual mutation probabilities q i . Second, for any pair of leaves u, v , the PHS distribution does not depend on the full structure of the tree, but only on a single sufficient statistic: the birth time of their LCA. This simplifies the computation of our proposed test statistic (described in the next section), and facilitates its theoretical analysis (see Section 6 ). Algorithm 1: Construction of the PHS test statistic Download figure Open in new tab 4.2 A PHS-based test statistic We are now ready to present our proposed PHS test statistic for the hypothesis problem in (3). For simplicity, in this section we assume the model parameters λ and q , as well as the birth times of all inner nodes of 𝒯, are known. As described in Appendix A.1, in our simulations we estimate λ and the birth times from the data. The collision probability q is assumed to be known also in simulations since, as detailed in Appendix A.1, q can be accurately estimated in various experimental settings. Moreover, as illustrated in Figure 10 (Appendix K), our test statistic is not sensitive to the exact value of q , and taking q = 1 /m , which corresponds to a uniform mutation distribution q i = 1 /m , works well. The calculation of our test statistic is outlined in Alg. 1. Given a candidate tree 𝒯 𝒮, it consists of two steps. Step 1 (calculating tail probabilities) . Based on Lemma 1 , for each pair of leaves ( u, v ) with latest common ancestor w = LCA( u, v ) in the reconstructed tree 𝒯, we compute the tail probability of its PHS value as if this tree were the ground-truth one, Each score s ( u, v ) can be viewed as a p-value (i.e., right-tail probability under ℋ 0 ) for the consistency of the observed sequences at the nodes u and v with the reconstructed tree. Since there are n leaves, Eq. (11) provides us with scores. The next step is to fuse these multiple p-values into a single score. Step 2 (fusing the multiple scores) . Motivated by multiple hypothesis testing procedures, we sort the N scores in increasing order, s 1 ≤ s 2 ≤ … ≤ s N , and focus on the first few smallest ones. Similar to the Benjamini-Hochberg (BH) procedure [ BH95 ] (see Remark 3 below for more details), for each i ≤ N we compute an adjusted score Finally, our test statistic is the minimal adjusted score, Note that cPHS(𝒯 𝒮) ≤ 1, since . Given a reconstructed tree 𝒯 𝒮, we reject ℋ 0 if cPHS(𝒯 𝒮) is below some threshold t < 1, and accept it otherwise. Remark 2. Let us provide some motivation and context for the above procedure. Often, some parts of the reconstructed tree are close to the ground truth. By construction, the corresponding scores s ( u, v ) in those parts are distributed approximately uniformly over [0, 1]. In contrast, other parts of the reconstructed tree may be less accurate. In those cases, some of the scores s ( u, v ) are expected to be closer to zero. Typically, the number of s ( u, v ) values that are close to zero is relatively small w.r.t. N ≫ 1. Hence, our setting is similar to high-dimensional multiple hypothesis testing problems, whereby only a few out of many hypotheses deviate from the null. Remark 3 (Relation to the BH procedure). As described in [ BHY09 , Section 2(b)], the BH procedure can be presented in terms of adjusted p-values, defined as Under this formulation, all hypotheses i for which are rejected. In our case, we do not aim to reject or accept individual hypotheses, corresponding to individual p-values . Instead, we adopt a more stringent criterion: if even a single adjusted score falls below the threshold, the entire tree is rejected. Accordingly, our test statistic coincides with the minimum BH-adjusted p-value: As seen in Figure 1(d) and further illustrated in Section 5 and in Appendix K, cPHS is highly powerful in distinguishing between accurate and inaccurate trees. Furthermore, as cPHS separates accurate and inaccurate reconstructed trees by a large margin (see Figure 1 ), it does not require a delicate tuning of the threshold t . In practice, as illustrated in Section 5 , a threshold value that depends only on the cutoff value ϵ of Eq. (3) works well for a wide range of model parameters (such as the sequence length k , the number of mutation states m , etc.). Remark 4 (Non-uniform mutation rate). Lemma 1 , and thus Eq. (11) , assume that the mutation rate λ is the same for all k characters. However, as long as characters are assumed to evolve independently, our PHS approach can be easily extended to the case of non-uniform mutation rates. Complexity and runtime Given a reconstructed tree 𝒯 𝒮, computing its cPHS value can be done in a polynomial number of operations. Specifically, as the algorithm passes over all leaf pairs for each character, its time complexity is 𝒪 ( kn 2 ). With our Python implementation, calculating cPHS for a tree with n = 1000 leaves and sequences of length k = 50 takes approximately two minutes on a standard PC. Missing data In practical settings, some entries of 𝒮 n may be missing, due to either limited sensitivity (the so-called stochastic drop out ) or heritable mechanisms (e.g. concomitant resection of two adjacent target sites [ JKQ + 20 ]). As described in Appendix B, our PHS approach can be easily extended to handle missing data, and empirically, it continues to perform well also in such cases. 5 Simulations We demonstrate through a comprehensive set of simulations the capability of the cPHS test statistic (13) to resolve the hypothesis testing problem outlined in (3). We compare cPHS with the three test statistics introduced in Section 3 : parsimony, parsimony tail-probability, and likelihood. In this section, the matrix 𝒮 n is fully observed; simulations with missing data appear in Appendix B. The simulations conducted in this section mirror those described in Section 3 , with each simulation including 1000 random ground-truth trees 𝒯 GT 𝒮 GT ; see also Appendix J. In all figures, 95% confidence intervals were computed via bootstrap. The observed data n from each ground-truth tree was given as input to seven reconstruction algorithms: Neighbor-Joining [ SN87 ], Cassiopeia-Greedy [ JKQ + 20 ], Shared-Mutation-Joining [ WZKY23 ], MaxCut and MaxCut-Greedy [ SR06 , JKQ + 20 ], Spectral and Spectral-Greedy solvers [ JKQ + 20 ]. The seven reconstructed trees were then classified as either accurate or inaccurate according to (3), using normalized RF ( d = d RF ) and parsimony ( d = d P ) as distance functions. Results for the other two distance functions, triplets ( d = d tri ) and likelihood ( d = d L ), are provided in Appendix K. Given a distance function and a cutoff value ϵ , for any tree denote by y ∈ {accurate, inaccurate} its true label according to (3). Similarly, denote by ŷ t the label estimated by some test statistic given a threshold t . For example, for cPHS, We first examine the performance of the four test statistics across all possible thresholds, then focus on specific choices of t . In the first simulation, we compute the true positive rate (TPR) and false positive rate (FPR) of each test statistic as a function of the threshold t . These rates are defined as follows: Plotting the TPR and FPR for all threshold values t gives ROC curves, presented in Fig. 2 . Panel (a) displays the results for d = d RF with ϵ = 1 / 4, and panel (b) displays the results for d = d P with ϵ = 0.01. It is apparent that, regardless of the selected threshold t , the performance of cPHS is superior to that of the other test statistics. Notably, the ROC curve of cPHS exhibits a steep ascent on the left side, indicating that cPHS achieves a satisfactory TPR even at low FPR levels. Download figure Open in new tab Figure 2: ROC (top) and Precision-Recall (bottom) curves of several test statistics for the hypothesis testing in (3). The distance measure is normalized RF with ϵ = 1 / 2 (left) and parsimony distance with ϵ = 0.01 (right). Dashed lines represent the expected performance of random guessing. The rates are calculated over 1000 realizations. Panels (c) and (d) of Fig. 2 present Precision-Recall curves. These metrics are defined as follows: A clear distinction in performance between cPHS and the other test statistics is evident under these metrics as well. Specifically, when the hypotheses in (3) involve the parsimony ( d = d P ), as seen in panel (d), cPHS achieves perfect Precision at a Recall value as high as 0.7. The simulations above were performed with sequence length k = 50 and number of mutated states m = 50. Next, we examined the performance as these parameters are varied. Figure 3 shows the AUC (area under the ROC curve) as a function of the sequence length k (left panel) and of the number of mutated states m (right panel). Additional results for varying values of n and ρ (bijectively related to λ , see (1)) can be found in Appendix K. As shown, under a wide range of settings, cPHS consistently outperforms the other test statistics by a significant margin. Download figure Open in new tab Figure 3: AUC values for the ROC curves of several test statistics for the hypothesis testing in (3) with d = d RF and ϵ = 1 / 2, as a function of the sequence length k (left), and as a function of the number of mutated states m (right). Each point is based on 1000 realizations of the simulation. We now turn to results obtained for specific choices of the threshold t . Our goal is to evaluate the robustness of cPHS when using a fixed threshold, rather than optimizing t for each parameter setting. To this end, we focus on the balanced accuracy metric, which is commonly used to evaluate classification performance. It is defined as the average of sensitivity and specificity: For each parameter configuration, we compute the balanced accuracy across several thresholds and record the maximum value. We then compare this optimal performance to that of cPHS using a fixed threshold t chosen based only on ϵ , independent of model parameters like k or m . Figure 4 summarizes this comparison. Specifically, the left panel shows results for different values of k , and the right panel for different values of m . In both panels the cutoff value is fixed at ϵ = 1 / 2, and the threshold of cPHS is fixed at t = 10 −3 . Two conclusions can be drawn from these results: First, the performance of cPHS is similar for a fixed threshold (based solely on ϵ ) and an optimally tuned one. Second, whether using fixed or optimal thresholds, cPHS consistently outperforms all other test statistics, evaluated with their own optimal thresholds. Additional results in Appendix K confirm that cPHS maintains strong performance with a fixed threshold across a wide range of values of n and ρ , as well as for a different value of ϵ . Download figure Open in new tab Figure 4: Balanced accuracy of several test statistics for the hypothesis testing in (3), with the same setting as in Figure 3 . Notably, cPHS with a fixed threshold achieves a significantly higher balanced accuracy than other approaches with optimally tuned thresholds. Finally, we illustrate the statistical power of cPHS over parsimony and likelihood by a different type of analysis. Specifically, we consider the following question: Suppose we are given a set of candidate trees, which includes trees reconstructed by various algorithms as well as the ground-truth tree. Is it possible to detect which tree is the ground-truth one? In other words, what is the probability that the test statistic of the ground-truth tree is better than its values for all reconstructed ones? Figure 5 illustrates results for various test statistics as a function of the number of mutated states m and of the mutation probability at a leaf ρ . Notably, the ground-truth tree is very rarely the most parsimonious or the maximum likelihood tree. In contrast, in approximately 40% of the cases (and up to 55% at high mutation rates), cPHS(𝒯 GT 𝒮 GT ) > cPHS(𝒯 𝒮) for all reconstructed trees 𝒯 𝒮. These results are in agreement with those described in Section 3 . Download figure Open in new tab Figure 5: Ability to detect the ground-truth tree as a function of the number of mutated states m (left) and of the mutation probability at a leaf ρ (right), with the same setting as in Figure 3 at k = 50. Each point is calculated over 1000 realizations. 6 Theoretical Support for PHS In this section, we present theoretical support for our PHS-based approach. In Section 4.2 we proposed to assess the accuracy of a full tree by the cPHS test statistic (13). This quantity, however, is difficult to analyze theoretically, as it depends in a non-trivial way on the tail probabilities of the PHS values at all leaf pairs in the tree. Instead, we analyze a simpler PHS-based test statistic, tPHS(𝒯 𝒮), defined below. Our main result is that even this simpler measure is significantly more powerful than parsimony in detecting inaccurate tree topologies. The PHS-based test statistic we consider in this section is defined as follows. Definition 4 (Total PHS). The total PHS of a full tree 𝒯 𝒮 is the sum of the PHS values over all its pairs of terminal leaves, where phs( u, v ) is defined in Eq. (9) . To simplify the analysis, we assume that the ground-truth tree has a homogeneous (binary) tree topology. Definition 5 (Homogeneous tree topology). A 𝒯 tree is called homogeneous if all its edges have the same length (elapsed time). The depth d of the tree is the number of edges from the root to any of its leaves. We consider the following scenario. Given the sequences at the terminal leaves, we assume that an algorithm reconstructed a tree topology 𝒯 which is nearly identical to the ground-truth 𝒯 GT , and differs from it by a single swap of two leaves u, v from different sides of the original tree. Our key result is that tPHS is qualitatively more powerful than parsimony in detecting such incorrect topologies. We show that the difference in tPHS between the incorrect tree and the ground-truth one, normalized by the tPHS standard deviation, is much larger than the parsimony difference, normalized by its standard deviation. In simple words, the tPHS measure can provably distinguish between such correct and incorrect reconstructed trees, whereas parsimony cannot. Next, we formally define a leaf swap. An illustration appears in Fig. 7 of Appendix F. Definition 6 (Leaf swap). Let u, v be a pair of leaves in a full tree 𝒯 GT 𝒮 GT . Denote by (𝒯 𝒮) u ↔ v the tree obtained from 𝒯 GT 𝒮 GT by swapping the leaves u and v . This swap operation is followed by a correction of the sequences in the ancestors of u and v to satisfy the non-modifiability constraint: at each character i ∈ [ k ], if a mutated ancestor w of u or v has an unmutated descendant, or if two of its descendants are mutated to different states, then the ancestor is set to be unmutated, . Further, denote the difference in tPHS and parsimony values following the swap by The following theorem shows that the difference is much larger than , both properly normalized.The proof appears in Appendix H. Theorem 1. Let 𝒯 GT 𝒮 GT be a homogeneous full tree of depth d whose sequences, of length k, were generated according to the model in Section 2.1 , with mutation rate λ, collision probability q, and mutation probability at a leaf ρ given in (1). Let u and v be a pair of leaves whose LCA ( u, v ) is the tree root. Let and be the respective change in tPHS and in parsimony, as defined ] in ?(15), following the swap u ↔ v. Then, for a sufficiently large d, their normalized expectations and satisfy the following lower bound and upper bound, respectively , The theorem holds provided that d is sufficiently large, with the required bound depending on λ . For instance, when λ ≤ 3, it suffices to have d ≥ 10. Let us discuss the consequences of this theorem. Consider a large tree with n ≫ 1 leaves, and depth d = log 2 n . Viewing the model parameters λ, q as fixed, according to Eq. (16) , for the change in tPHS to be significant it suffices That . That is, the required sequence length scales only polylogarithmically with the number of leaves, k ≳ (log 2 n ) 2 . We note that, as proven in [ ESSW99 ] and [ WZKY23 ], a similar relation between k and n suffices for accurate tree reconstruction under reversible and non-modifiable mutation models, respectively. In contrast, for the change in parsimony to be significant, a necessary condition is that the right-hand side of (17) is large. This translates into a significantly longer sequence of length k ≳ n 2 / (log 2 n ) 5 , which grows polynomially with n . In simple words, compared to tPHS, parsimony requires exponentially longer (in tree depth) sequences to detect an incorrect topology of the form of a leaf swap. Remark 5 (leaf swap on the same side of tree). Theorem 1 addresses the swap of leaves whose LCA is the tree root. The theorem can be easily generalized to the swap of leaves u, v whose LCA( u, v ) is an inner node. This would result in two changes to the RHS of Equations (16) and (17) . First, instead of d we would have d − l , where l is the depth of LCA( u, v ). Second, instead of k we would have k · e − λl/d , the expected number of unmutated characters at the node LCA( u, v ). As long as l is much smaller than d , namely LCA( u, v ) is close enough to the root, the consequences of the theorem continue to hold. Specifically, tPHS would be able to detect such a leaf swap with a sequence length polylogarithmic in n , whereas parsimony would require a sequence length polynomial in n . Remark 6 (several leaf swaps). Theorem 1 considers a single leaf swap. The theorem can be readily extended to the case of several leaf swaps, provided that the different LCAs of the leaf pairs are located at distinct branches of the tree. In particular, the LCAs must be inner nodes; see Remark 5. In this case, the effect of the leaf swaps on both the tPHS and parsimony is additive. 7 Discussion In this work, we introduced a homoplasy-based approach to distinguish between accurate and inaccurate reconstructed trees. As demonstrated via simulations in Section 5 and supported by theoretical analysis in Section 6 , our method offers significantly greater discriminative power than traditional parsimony. Intuitively, this advantage may be explained as follows: parsimony aggregates a local measure — the number of mutations per edge — without accounting for the broader tree structure. In contrast, PHS integrates information across all pairs of leaves, both close and distant, thereby capturing the consistency of the tree’s topology as a whole. In our work, we considered the minimal adjusted p-value as our test statistic (13). The key parameter of our procedure is the threshold t against which the cPHS statistic is compared. In practical settings, a suitable threshold t can be estimated via simulations designed to reflect the characteristics of the experimental setting. Specifically, one may generate multiple ground-truth trees, for example via Cassiopeia [ JKQ + 20 ]. Given the observed sequences for each ground-truth tree, various reconstruction algorithms can be run, followed by computing the cPHS scores for their outputs. A suitable threshold can then be determined based on the cPHS scores obtained for the most accurate reconstructions. We remark that there are other possible schemes to fuse the scores in (11) into a single test statistic, e.g. the harmonic mean [ Wil19 ]. Our simulations indicate that the statistical power of the harmonic mean is comparable to that of cPHS. However, it is less practical, as in contrast to cPHS, the harmonic mean is very sensitive to the specific value of the threshold t . In this manuscript, we focused on a non-modifiable model of evolution, whereby Lemma 1 , and consequently our cPHS test statistic (13), are based on this assumption. Our proposed homoplasy approach may also be beneficial for assessing the accuracy of reconstructed trees under reversible models. This requires deriving an analogue of Lemma 1 for a given reversible model of evolution. We leave this extension for future work. Another interesting direction for future research is to use a homoplasy approach not merely to assess the accuracy of trees, but in fact to develop a PHS-based tree reconstruction algorithm. As we illustrated both empirically and theoretically, PHS provides a better tree accuracy measure than parsimony. Hence, it would be interesting to explore if a PHS-based algorithm can reconstruct more accurate trees than those found by current methods that aim to find the tree with maximum parsimony for the observed data. A Handling Model Uncertainty and Incomplete Reconstructions In this section, we address some practical aspects in applying the cPHS statistic. Specifically, we discuss how to estimate unknown model parameters required for the computation; how to handle cases where only the reconstructed tree topology is available without internal node sequences; and how to prevent saturation of the cPHS score in reconstructed trees with very few homoplasies. A.1 Unknown model parameters and birth times As discussed in the main text, our PHS test statistic depends on two model parameters: the mutation rate λ and the collision probability q . In addition, in terms of the given tree, it depends on the birth times of the inner nodes. In this section, we discuss their estimation. Let be the fraction of mutated characters in the sequences at the observed cells. Then, as in [ JKQ + 20 ], the mutation rate λ can be accurately estimated from via the plug-in estimate ; see (1). The collision probability q , in contrast, cannot typically be accurately estimated given the result of a single or few clones, since we do not know when up in the tree an observed mutation occurred. However, this parameter is constant for a fixed experimental system [ JKQ + 20 , SS22 ]. It is thus possible to estimate it from single-cell RNA-seq data of many clones or even from bulk RNA-sequencing of the target sites, as done in [ JKQ + 20 , SS22 ]. In our simulations, we assume knowledge of q , and estimate λ from the observations. Finally, the birth times τ of the inner nodes are estimated using the maximum-likelihood-based branch length estimator proposed in [ PRYS23 ]. A.2 Reconstructed topology 𝒯 without internal node sequences 𝒮 In the main text, it was assumed that the output of a given algorithm is a reconstructed full tree 𝒯 𝒮, which includes both the tree topology 𝒯 and the internal node sequences 𝒮. However, some algorithms output only a reconstructed tree topology without a corresponding. Given 𝒯 and 𝒮 obs , there is more than one possible set of internal node sequences 𝒮 consistent with the non-modifiability model assumption. Since to compute the cPHS score requires also the set of internal sequences, the cPHS score “induced” by only a reconstructed topology is not uniquely defined. Similarly, the hypothesis testing (3) with either the parsimony or the likelihood distance is also not well-defined. In this section, we propose a procedure to reject candidate topologies w.r.t. the parsimony distance , even without internal sequences. We first present the procedure, and then its justification. Assuming we are given only a candidate topology 𝒯 without reconstructed internal sequences, we first compute the solution to the small MP problem 𝒮 MP , under the non-modifiability constraint. This is described in Appendix C and Alg. 2. Next, we compute the test statistic cPHS (𝒯 𝒮 MP ). If our test rejects 𝒯 𝒮 MP , then we reject the topology 𝒯 itself. The justification of this procedure is based on the following lemma, whose proof appears in Appendix E. It states that the solution to the small MP problem maximizes the cPHS score among all the possible sets of sequences that comply with the non-modifiability constraint and coincide with the observed sequences at the leaves. Lemma 2. Let 𝒯 be a fixed tree topology, and let 𝒮 n be a given set of sequences at its leaves with no missing data. Denote by 𝒮 MP a solution of the small MP problem. Then for any set of sequences 𝒮 at all tree nodes that complies with the non-modifiability constraint and whose leaves coincide with 𝒮 n , it holds that Suppose that our procedure rejected 𝒯 𝒮 MP . Assuming our test result is correct, 𝒯 𝒮 MP satisfies hypothesis ℋ 1 of (3), namely 𝒯 𝒮 MP > (1 + ϵ ) · M (𝒯 GT 𝒮 GT ). Hence, by definition of the small MP problem, 𝒯 𝒮 MP satisfies ℋ 1 of (3) for any 𝒮. The topology 𝒯 should thus be rejected, as required. A.3 Very few homoplasies In an extreme case where the reconstructed full tree has very few homoplasies, it may occur that all sorted scores satisfy In this case, using (13) leads to a saturation in cPHS, outputting the maximal possible value of 1. This precludes a finer comparison of different trees that both obtain cPHS = 1. To overcome this issue, we exclude leaf pairs with no homoplasies, { i : s i = 1}, before taking the minimum. Formally, Eq. (13) is modified to B Missing Data In the main text, we assumed a complete (i.e., fully observed) data matrix 𝒮 n . In practice, however, some entries in 𝒮 n , corresponding to characters at the tree leaves, may be unknown. The CRISPR-Cas9 evolution model has two distinct mechanisms of missing data. The first is stochastic missingness, which arises from imperfections in the sequencing process, leading to randomly missing entries in 𝒮 n . The second is heritable missingness, which occurs when characters at internal nodes are absent due to transcriptional silencing or double Cas9 resection. These missing characters behave analogously to mutations: their unknown state is non-modifiable and is inherited by all descendant nodes, including the leaves. Let p sto and p her denote the rates of stochastic and heritable missingness in 𝒮 n , respectively. Assuming independence, the overall missing data rate is then given by p miss = p sto + p her p sto − · p her . To deal with the presence of missing data, two minor modifications are required in the definition of the cPHS measure. Let u and v be a pair of leaf nodes. First, we revise the definition of phs i ( u, v ) in Eq. (8) as follows: if the i -th character of either u or v is missing, then phs i ( u, v ) = 0. Second, we modify Eq. (11) in Step 1 of the cPHS computation procedure by replacing k with k eff ( u, v ), defined as the number of characters that are non-missing in both u and v . As illustrated in Figure 6 , the adapted cPHS performs well in the presence of missing data, and in particular, it retains its advantage over the other test statistics. Interestingly, the performance of cPHS, as well as that of the other test statistics, improves as p miss increases. This trend, however, may be attributed to a technical factor: the balance between accurately and inaccurately reconstructed trees shifts with increasing p miss . Specifically, as p miss increases, the reconstruction task becomes more challenging, leading to a higher proportion of inaccurate trees produced by the algorithms. It is reasonable to assume that this shift is the primary driver of the observed improvement in the test statistics’ performance. Download figure Open in new tab Figure 6: AUC values (left) and balanced accuracy (right) of several test statistics for the hypothesis testing in (3), with the same setting as in Figures 3 and 4 at k = m = 50. The x-axis is the total proportion of missing data p miss . C Small and Large Parsimony Problems Let us recall the definitions of the small and large maximum parsimony (MP) problems [ Fel04 , Chapters 2 and 4]. Definition 7 (Small MP). In the small MP problem, given a tree topology 𝒯 and observed sequences at the leaves 𝒮 n , the goal is to reconstruct the sequences at the ancestral (internal) nodes such that the resulting tree has the smallest number of mutations: Algorithm 2: Solution to the small MP problem under non-modifiable model Download figure Open in new tab Whereas small MP optimizes only over the internal sequences for a fixed topology, large MP additionally optimizes over the tree topologies: Definition 8 (Large MP). In the large MP problem, given the observed sequences at the leaves 𝒮 n , the goal is to reconstruct a full tree with the smallest number of mutations: The small MP problem can be solved efficiently by the Sankoff algorithm [ San75 ]. The non-modifiability constraint can be forced by setting the substitution matrix accordingly. The resulting algorithm, presented in Alg. 2, is simple, and bears a resemblance to Fitch algorithm [ Fit77 ]. In the absence of missing data, the solution to the small MP problem under the non-modifiability constraint is unique [ SSCR23 ]. The large MP problem, in contrast, is NP-hard in general [ Fel04 ], and remains so under the non-modifiable model [ DJS86 , SSCR23 ]. D Parsimony Tail Probability In this section, we derive the parsimony distribution conditional on a given candidate tree topology, where the inner sequences are random. The distribution is calculated recursively along the tree topology Let w be an unmutated node in the tree, and denote the number of mutations at the i -th character in the subtree whose root is w by . If w is a leaf, then no mutation can occur from w downwards, namely . Next, consider the case where w is an inner node. Denote its children by u and v , and the edge length from w to u and to v by τ w,u and τ w,v , respectively. Recall that the edge length is the time elapsed between the birth time of a node and its child. Let U and V be the events of mutation occurrence at the i -th character from w to u and to v , respectively. Then and . In addition, by the Markovian property of the generative process, U and V are independent given that w is unmutated. By the law of total probability, If both u and v are un mutated (namely under . As and are independent, . If v is mutated but not due to the non-modifiability of the evolution process. Hence, . Similarly, . Finally, if both u and v are mutated, due to non-modifiability, and thus . Putting everything together, we conclude This recursive equation allows us to compute the parsimony distribution at a node w at the i -th character, , as a function of its children distributions, and . The distribution of the total parsimony (sum of mutations at all characters) of a tree topology 𝒯 with random sequences 𝒮 is given by convolving the individual distributions , The second equality follows from the assumption of uniform mutation rate; see Remark 4. Finally, the parsimony tail probability for a given tree topology 𝒯 is ℙ[ M (𝒯 𝒮) ≤ M root ]. E Proofs of Lemmas 1 and 2 Proof of Lemma 1 . Let p u,v = ℙ[phs 1 ( u, v ) = 1], where the probability is over the random sequences 𝒮 GT . Then, by the definition of homoplasy, with w = LCA( u, v ), Let us analyze each of the three terms in the RHS of (22). The first expression is simply the probability that the first character in s ( w ) is unmutated, which is given by Regarding the second term, since characters at node u and at node v evolve independently, conditional on their LCA w being unmutated, The probability that the character at u is unmutated given that w is unmutated is where τ w,u = τ u − τ w = 1 − τ w , where the second equality follows from the fact that u is a leaf. Hence, . As a similar argument holds for v , As for the third term, since a mutated character is i with probability q i , the probability of a collision (same character in both nodes) is given by Inserting all of these into (22) yields . The Binomial distribution in (10) of the lemma follows from the fact that characters at different locations evolve independently with the same parameters λ and q . Proof of Lemma 2 . Let 𝒮 be a set of sequences that coincides with 𝒮 GT at the leaves. Assume that for any pair of leaves u, v , Then, by (11), s ( u, v |𝒮 MP ) ≥ s ( u, v | 𝒮). Therefore, upon computing (10), (11), sorting the s -values, adjusting them by (12) and taking the minimum as in (13), (18) of the lemma follows. Hence, we now prove that (23) indeed holds. Let u, v be a pair of leaves, and let w = LCA( u, v ). For simplicity, we focus on the first character location. Since characters evolve identically and independently, the same argument holds for all the characters. Let z 1 , z 2 be the immediate children of w . There are four generic possible scenarios for and (1, 1). All the other possibilities are equivalent to one of these four. Under the non-modifiable model, must be 0 in the first three cases. The only freedom of choice left to the inner sequences reconstruction algorithm is in the case , which we now analyze in detail. In this case, all the descendants of w have 1 in their first character. Hence, assigning minimizes the PHS of all leaf pairs that are descendants of w , including u and v . It is left to show that in 𝒮 MP , w is indeed assigned with 1 in this case. Let us analyze the effect of on the total number of mutations, M (𝒯 𝒮 MP ). Assigning increases the total number of mutations by either 0 (if the sibling of w has 1) or 1 (if it has 0). Assigning w with 0, on the other hand, increases the total number of mutations by at least 2, as two changes are required from to and to . Hence, the small MP solution assigns with 1, as required. F Illustration of a Leaf Swap, and its Effects on PHS and Parsimony Let us illustrate the difference in tPHS and parsimony values following a leaf swap. For simplicity, the sequence length in our example consists of a single character, k = 1. Suppose that in the ground-truth tree, the immediate offspring of the root mutated to two different states: the left branch mutated to 1 and the right branch to 2. As a result, the state of all the leaves in the left subtree is 1, and it is 2 in the right subtree; see Figure 7 (left) for an illustration. Now, suppose that a tree reconstruction algorithm made a single error of swapping the leftmost and rightmost leaves; see Figure 7 (right). In this example, there are no homoplasies in the ground-truth tree, and all its leaf pairs have phs( u, v ) = 0. In the reconstructed tree, in contrast, a non-negligible fraction of leaf pairs have a PHS of 1. As stated in Lemma 10 below, in this case the increase in tPHS scales as the number of leaf pairs, which is roughly 4 d . This is comparable to the tPHS standard deviation over random inner sequences; see Lemma 4 . Hence, PHS is sufficiently powerful to identify the reconstructed tree as incorrect. In contrast, as presented in Lemma 5 , the parsimony increase is only linear in d . This is negligible compared to the parsimony standard deviation over random inner sequences, which scales as the number of leaves, n = 2 d ; see Lemma 3 . As a result, the parsimony measure is unable to distinguish between the correct and the incorrect tree. Download figure Open in new tab Figure 7: Effect of a leaf swap. Left: a ground-truth tree where mutations occurred early in the tree (at the green nodes). As a result, each of its two subtrees has a different mutated state. Right: a reconstructed tree which differs from the GT tree by a single swap (the swapped leaves are marked in orange). The inner node sequences are reconstructed to satisfy the non-modifiability assumption. G Auxiliary Lemmas for Theorem 1 The proof of Theorem 1 makes use of four auxiliary lemmas, outlined below. In the following, ℒ (𝒯) denotes the set of leaves of a tree 𝒯. G.1 Mean and variance of parsimony and tPHS We start with two lemmas, regarding the mean and variance of parsimony and tPHS. The first lemma provides exact expressions for the parsimony mean and variance, assuming a homogeneous tree. The second lemma provides upper bounds for the mean and variance of tPHS in general trees, not necessarily homogeneous. Their proofs appear in Sections I.1 and I.2 , respectively. Lemma 3. Let 𝒯 GT be a homogeneous tree of depth d, whose root node is unmutated. Let 𝒟 be the distribution of its set of sequences 𝒮 GT of length k, generated according to the non-modifiable model of Section 2.1 with mutation rate λ. Then, taking expectation over 𝒟, for α d = e − λ/d ≠1 / 2, and Lemma 4. Consider the same setting as in Lemma 3 , but with an arbitrary tree 𝒯 GT (not necessarily homogeneous) . Denote by its number of leaf pairs. Then, taking expectation over 𝒟, and where q is the collision probability and ρ is the probability of observing a mutation at a leaf given in (1). Comparing Lemmas 3 and 4 reveals an important difference between parsimony and tPHS: the parsimony does not depend on the collision probability q . In contrast, the mean and variance of tPHS are linear in q . As a result, in general, a small value of q - which is directly related to the number of homoplasies - increases the statistical power of cPHS; see (16). Remark 7. Lemma 3 assumes that α d ≠1 / 2. In the case α d = 1 / 2, the mean and variance read This can be easily verified by plugging α d = 1 / 2 in the proof of the lemma, or by taking the limit α d → 1 / 2 in Eqs. (24) and (25) . G.2 The effect of a leaf swap on parsimony and on tPHS We present two lemmas regarding the mean change in parsimony and tPHS, following a swap of two leaves u, v from different sides of the tree (see Section 6 ). For convenience, the lemmas consider the change in parsimony and tPHS at a single character. Since both parsimony and tPHS are additive in the sequence length k , the mean change in these quantities is simply k times that of a single character. Definition 9. [per-character parsimony] Given a full tree 𝒯 𝒮, for any i ∈ [ k ] denote by M i (𝒯 𝒮) the number of mutations in the i -th character across all edges of the tree. Definition 10. [per-character tPHS] Given a full tree 𝒯 𝒮, denote by tPHS i (𝒯 𝒮) the number of leaf pairs with a PHS of 1 at their i -th character, i ∈ [ k ]: The factor of 1 / 2 compensates for double summation of pairs. Note that by definition, the parsimony and Tphs of a tree are the sum of its per-character parsimony and tPHS scores: and , respectively. In addition, tPHS i (𝒯 𝒮) = k · ∑ ( u,v ) phs i ( u, v ). Let us now present the lemmas. Their proofs appear in Appendices I.3 and I.4, respectively. Lemma 5. Consider the same setting as in Lemma 3 . Further suppose λ ≤ d/ 2. Let u, v ∈ ℒ (𝒯 GT ) be two leaves whose LCA is the tree root. Then, taking expectation over 𝒟, for any i ∈ [ k ], the expected parsimony change in the i-th character following the leaf swap u ↔ v is given by Lemma 6. Assume the conditions of Lemma 5 . Then, taking expectation over 𝒟, the expected tPHS change in the i-th character following the leaf swap u ↔ v satisfies Lemmas 5 and 6 demonstrate a key difference between parsimony and tPHS: Following the leaf swap described above, the mean tPHS value increases exponentially with respect to the tree depth d , whereas on average, parsimony increases only linearly. H Proof of Theorem 1 Proof . First, we prove the lower bound (16) on the normalized change in tPHS due to a leaf swap. Let α d = e − λ/d . For any d ≥ 2 λ , Inserting this inequality into (31) of Lemma 6 yields Next, we apply Lemma 4 . Since the tree is homogeneous, its number of leaves is n = 2 d . Now, for . Hence, Eq. (27) of the lemma gives Equation (16) follows by combining this with (33). Next, we prove the upper bound (17) regarding the normalized change in parsimony. Combining the additivity of the parsimony in the sequence length k with (31) of Lemma 5 regarding the change in parsimony at a single character, gives that where the inequality above follows from 1 − α d ≤ λ/d . Next, we apply Lemma 3 . For a sufficiently large d, α d > 1 / 2, and thus satisfies the assumption of the lemma. Further, for a sufficiently large d , the expression in the curly brackets in the RHS of (25) is larger than 2 / 3. Hence, where the second inequality above follows from α d ≤ 1. Combining this with Eqs. (1) and (32) yields Equation (17) follows by combining this with (34). I Proofs of Lemmas 3 , 4 , 5 and 6 To prove the lemmas, we first introduce some definitions. Recall that ℒ (𝒯) denotes the set of leaves of a tree 𝒯. For convenience, for any leaf z define phs( z, z ) = 0. In addition, we define the following three quantities regarding the origin of observed mutations, that will be used in our proof. See Fig. 8 for their illustration. Download figure Open in new tab Figure 8: Illustration of the origin of an observed mutation (Def. 11) at a node u (colored in orange) with k = 1. The most distant node with the same mutation, A 1 ( u ), is colored in green node. Their edge distance is E 1 ( u ) = 2. Definition 11. Let 𝒯 𝒮 be a full tree whose sequences were generated by the non-modifiable mutation process of Section 2.1 . Let u ∈ ℒ (𝒯) be a leaf in the tree. For each character location i , we denote by A i ( u ) the most distant ancestor of u that had a mutation in the i -th character. If either or but its father node was unmutated, then we define A i ( u ) = u . In addition, we denote by E i ( u ) the number of edges from A i ( u ) to u , and by the subtree whose root is A i ( u ). Note that by non-modifiability of the mutation model, . I.1 Proof of Lemma 3 (parsimony mean and variance) The proof of Lemma 3 makes use of the following auxiliary lemma, to be proved shortly. Lemma 7. Assume the same setting as in Lemma 3 . Let 𝒯 l be a subtree of 𝒯 GT of depth l ≤ d. Let D l be the distribution of its set of sequences 𝒮 l ⊆ 𝒮 GT , conditioned on the event that the sequence at the root node of 𝒯 l is unmutated. For any i ∈ [ k ], denote the parsimony in the i-th character of the subtree 𝒯 l 𝒮 l by . Then, taking expectation over 𝒟 l , and Proof of Lemma 3 . By the linearity of expectation, , where is defined in Lemma 7 . As characters at different locations i ∈ [ k ] evolve independently from each other, The lemma follows by taking l = d in Eqs. (35) and (36) of Lemma 7 . Proof of Lemma 7 . Without loss of generality, we may assume that i = 1. To simplify notation, we omit the subscript 1 from . We prove the claim by induction on l . For l = 0, the subtree consists of only the root, and thus M 0 = 0 deterministically. As the RHS of Eqs. (35) and (36) vanish at l = 0, the induction base is proved. Next, assume that the induction hypothesis holds, namely (35) and (36) hold for some l − 1 ≥ 0. We now prove it holds for l . Let E l = 𝔼[ M l ] and V l = 𝕍[ M l ] be the mean and variance, respectively, of the number of mutations in a subtree of depth l whose root has an unmutated sequence. The mean and variance depend only on l and not on the specific subtree, because all subtrees of a homogeneous tree are homogeneous as well. Let and be the parsimony scores of the subtrees at the left and right branches of 𝒯 l 𝒮 l , respectively, each conditioned on the event that the corresponding subtree roots are unmutated. For future use, note that by the Markovian property of the generative process, and are independent given that the root is unmutated. To analyze E l and V l , we split into cases. Let A ∈ {0, 1, 2} be the number of mutations from the root of 𝒯 l to its immediate children; see Figure 9 for an illustration of the different possible events. We first calculate the mean and variance conditioned on a specific value of A . Under the event A = 0 (case (i) in Figure 9 ), Download figure Open in new tab Figure 9: Given an unmutated root, there are four possible scenarios for the states of its two children: (i) no mutation in both branches, (ii) a mutation only in the left branch; (iii) a mutation only in the right branch, and (iv) mutations in both branches. Here X, Y ∈ [ m ]. and Conditioned on A = 1, there are two scenarios: the mutation can be either in the left or in the right branch (cases (ii) and (iii) in Fig. 9 , respectively). Due to the non-modifiability of the mutation process, no further mutations can occur in the mutated branch. As a result, in case of a mutation in the left branch, and in case of a mutation in the right branch. Since the scenarios of a mutation at the left and at the right branches have equal probabilities, and Conditioned on A = 2 (case (iv) in Figure 9 ), due to the non-modifiability of the mutation process, The probability of observing a mutation along a single edge, given that the parent node is unmutated, is 1 − α d . By independence of and , we have that ℙ[ A = 2] = (1 − α d ) 2 , ℙ[ A = 1] = 2 α d (1 − α d ), and . Hence, by the law of total expectation, Together with Eqs. (37a) , (38a) and (39) , we obtain Inserting the induction hypothesis (35) for l − 1 yields that (35) hold for l . Similarly, by the law of total variance, The second term above satisfies so that Inserting the expressions derived for the mean and variance conditioned on A , Eqs. (37), (38) and (39) give that where the last equality follows from (35). Plugging in the induction hypothesis (36) at l − 1 yields (36) at l . I.2 Proof of Lemma 4 (PHS mean and variance) Proof . Let u, v ∈ ℒ (𝒯 GT ) be a pair of leaves in the tree, and denote their LCA by w . Let τ w ∈ [0, 1) be the birth time of w . By Lemma 1 , k · phs( u, v ) follows a Binomial distribution with k trials and success probability . Hence, using (1), Equation (26) of the lemma follows from the definition of tPHS (14) and the linearity of expectation. Next, let us consider the variance. By (10), for a single pair of leaves, Equation (27) follows from the fact that for any collection of identically distributed random variables X i . I.3 Proofs of Lemma 5 (mean parsimony change following a leaf swap) For a fixed tree 𝒯 GT , the change in parsimony following a leaf swap between u, v is a random variable that depends only on the random quantities E i ( u ) and E i ( v ), defined in Def. 11. Indeed, as the following auxiliary lemma shows, conditional on E i ( u ) and E i ( v ), the change is deterministic. The proof of Lemma 5 follows by combining this result with Lemma 9 below, which characterizes the distribution of E i ( u ). The proofs of these auxiliary lemmas appear in Appendices I.5 and I.6. Lemma 8. Let 𝒯 GT 𝒮 ; GT be a full tree, and let u, v ∈ ℒ (𝒯 GT ) be two leaves whose LCA is the tree root. Then, for any i ∈ [ k ], the parsimony change in the i-th character following the leaf swap u ↔ v is given by Lemma 9. Let 𝒯 GT 𝒮 GT be a full tree with a homogeneous topology of depth d, generated according to the non-modifiable model of Section 2.1 with a mutation rate λ. Denote ρ = 1 − e − λ as in (1), and α d = e − λ/d . Let u be a leaf in the tree. For any i ∈ [ k ], if the i-th character of u is unmutated , , then E i ( u ) = 0. Otherwise, its distribution is given by We now prove the lemma. Proof of Lemma 5 . Denote and . Let Δ = M i ((𝒯 𝒮) u ↔ v )) − M i (𝒯 GT 𝒮 GT ). The quantity of interest is 𝔼[Δ]. To this end, denote the event A = { s u ≠ s v }. Under A c , Δ = 0. Hence, Next, we split A into three cases, depending on whether s u or s v were mutated, or both: By a symmetry argument, ℙ[ A 2 | A ] = ℙ[ A 3 | A ] and 𝔼[Δ | A 2 ] = 𝔼[Δ | A 3 ]. Hence, Next, we compute the probabilities in the equation above. Since the LCA of u and v is the root, which is unmutated by construction, s u and s v are independent random variables. We thus obtain Inserting this into (42) yields Next, we compute the conditional expectations. Denote e u = E i ( u ) and e v = E i ( v ). By combining this with the law of total expectation and (40) of Lemma 8 , Since s u = 0, E i ( u ) = 0 follows by definition. Hence, Combining this with (41) of Lemma 9 , we obtain By similar arguments, Let us calculate the sum in the RHS of (45), Plugging this into (45) gives that The lemma follows by combining (43) with (44) and (46). I.4 Proof of Lemma 6 (mean tPHS change following a leaf swap) Similar to parsimony, the change in tPHS following a leaf swap between u, v is a random variable that depends only on the two random quantities E i ( u ) and E i ( v ), defined in Def. 11. Indeed, as the following auxiliary lemma shows, conditional on E i ( u ) and E i ( v ), the change is deterministic. The proof of Lemma 6 follows by combining this result with Lemma 9 above, which characterizes the distribution of E i ( u ). The proof Lemma 10 appears in Appendix I.5. Lemma 10. Let 𝒯 GT 𝒮 GT be a homogeneous full tree, and let u, v ∈ ℒ (𝒯 GT ) be two leaves whose LCA is the tree root. Suppose that either or . Then, for any i ∈ [ k ], We now prove the lemma. Proof of Lemma 6 . Denote and By similar arguments as in the proof of Lemma 5 , (43) holds for Δ = tPHS i ((𝒯 𝒮) u ↔ v ) − tPHS i (𝒯 GT 𝒮 GT ) as well. That is, where A 1 = { s u ≠ 0 ∧ s v ≠ 0 ∧ s u ≠ s v } and A 2 = { s u = 0 ∧ s v ≠ 0}. Since u and v are from different sides of the tree, s u and s v are independent, as well as e u and e v . By combining this with the law of total expectation and (47) of Lemma 10 , Invoking Lemma 9 gives that Observe that By symmetry, . Inserting these expression into (49) yields Similarly, by combining the law of total expectation, (47) of Lemma 10 and (41) of Lemma 9 , Inserting the last two results into (48) gives (31) and completes the proof of the lemma. I.5 Proofs of lemmas 8 and 10 (deterministic parsimony and PHS change following a leaf swap) The proofs of the lemmas use the following two propositions regarding the quantities defined in Def. 11. Proposition 1. Let 𝒯 𝒮 be a full tree, and let u be a leaf in the tree. Let i ∈ [ k ]. Then for any , In addition, let w be the parent of A i ( u ). Then Proof . To simplify the notations, for any node z , denote . Let . There are two cases to consider: s u = 0 or s u ≠ 0. If s u = 0, then by non-modifiability all ancestors of u are unmutated. Hence, A i ( u ) = u , and . In this case, (50) trivially holds. Next, suppose s u ≠ 0. Since A i ( u ) is an ancestor of u , whose character is mutated, by the non-modifiability . Since any is either A i ( u ) or an offspring of A i ( u ), (50) follows by the non-modifiability. Finally, (51) follows from the definition of A i ( u ) as the most distant mutated ancestor of u . Proposition 2. Let 𝒯 GT 𝒮 GT be a full tree, and let u be a leaf in the tree. Let i ∈ [ k ], and denote . Let 𝒯 𝒮 be the tree after setting and correcting for non-modifiability (cf. Def. 6). Then, for all ancestors z of u that belong to it holds that in the original tree , whereas in the modified tree . The sequences at all the other nodes, except u, are the same in 𝒯 GT 𝒮 GT and 𝒯 𝒮. Proof . Let i ∈ [ k ]. If , then the lemma holds trivially. Hence, we assume that there is at least one ancestor of u in , namely where w is the parent of u . For convenience, for any node z denote . By Proposition 1 , s z = x for any in the original tree 𝒯 GT 𝒮 GT . In the modified tree 𝒯 𝒮, by definition, s u = y . Since , the sibling of u , denoted ū , satisfies s ū = x ≠ y . Hence, by non-modifiability, s w = 0 in the modified tree. As a result, all ancestors of u are corrected from x to 0. Since violations of non-modifiability can only occur from the modified leaf and upwards, the proof is complete. We now prove the lemmas. Proof of Lemma 8 . For any node z , denote by the state of the i -th character of z in 𝒯 GT 𝒮 GT and by its state in (𝒯 𝒮 u ↔ v ). Denote Δ = M i ((𝒯 𝒮) u ↔ v ) − M i (𝒯 GT 𝒮 GT ). If s u = s v , then (𝒯 𝒮) u ↔ v = 𝒯 GT 𝒮 GT , so Δ = 0 and (40) holds. Next, we analyze the case s u ≠ s v . Since u and v are from different sides of the tree and the root sequence is fixed with unmutated states, we can decompose the parsimony change as where Δ u and Δ v are the parsimony changes due to the modification of s u (while keeping s v unchanged) and of s v (while keeping s u unchanged), respectively. Let us begin by calculating Δ u . Denote a u = A i ( u ) and . We split into cases: e u > 0 and e u = 0. First, suppose that e u > 0. According to Proposition 2 , the only inner nodes whose sequences are modified following the swap are ancestors z of u that satisfy . Specifically, they are modified from s u to 0. Hence, for each such node z , there is a new mutation along the branch to the child which is not u or one of its ancestors. The number of these nodes z is E i ( u ). In addition, since , there is no longer a mutation from a u ’s parent to a u . Finally, we need to ask if there is a new mutation from w , the parent of u , to u . Since , if s v = 0 then there is no such mutation. Suppose now s v > 0. Since , and thus . Since , there is a new mutation from w to u . In total, we obtain that if e u > 0, then where is the indicator of the event { s v > 0}. Under the assumption e u > 0, we have . Hence, we can write Next, suppose that e u = 0. Then, by definition of e u , s w = 0. In addition, since , Proposition 2 implies that . Recall that . Let us split into three cases: (i) s u = s v , (ii) s u = 0 and s v > 0, and (iii) s u > 0 and s v = 0. In case (i), and Δ u = 0. In case (ii), and Δ u = 1. Finally, in case (iii), and Δ u = −1. It can be thus verified that (53) holds also in the case e u = 0. By similar arguments, one can show that where is either 0 or greater than 0. Inserting this together with (53) into (52) yields (40). In the following proof, we denote the PHS of a leaf pair before and after the swap by and , respectively. Proof of Lemma 10 . Let i ∈ [ k ]. For any node z , denote by the state of the i -th character of s ( z ) in 𝒯 GT 𝒮 GT , and by its state in (𝒯 𝒮) u ↔ v . Denote Δ = tPHS i ((𝒯 𝒮) u ↔ v ) − tPHS i (𝒯 GT 𝒮 GT ). Further denote x = s u and y = s v . If x = y , then (𝒯 𝒮) u ↔ v = 𝒯 GT 𝒮 GT , in which case the LHS of (47) vanishes. Since, by assumption, x = y = 0, E i ( u ) = E i ( v ) = 0, and the RHS of (47) vanishes as well. Next, we analyze the case x ≠ y . For convenience, we denote by the contribution of a specific pair of leaves at their i -th character to the total PHS change Δ, Let and be the sets of tree leaves whose LCA is A i ( u ) and A i ( v ), respectively. Further let W = U ∪ V and R = ℒ (𝒯 GT ) \ W . Then where the union is over disjoint sets. Since for for any , we can decompose the total PHS change as Let U * = U \ { u } and V * = V \ { v }. Then U × V = ( U * × V * ) ∪ ({ u } × V * ) ∪ ( U * × { v }) ∪ ({ u } × { v }), U × U = ( U * × U * ) ∪ ({ u } × U * ) ∪ ( U * × { u }) ∪ ({ u } × { u }), and similarly for V × V . Therefore, where we used the fact that by definition, Δ( u, u ) = Δ( v, v ) = 0. Let us analyze each of the seven terms in the RHS of (54). First, we show that the first two terms vanish, since To prove (55), let , and consider three possible cases for the leaf z : either z ∈ R, z ∈ U , or z ∈ V . If z ∈ R , then also . By Proposition 2 , the sequences at both , and are not affected by the leaf swap, and thus . Next, suppose that z ∈ U . By Proposition 1 , . By Proposition 2 , also , that is . Hence also in this case. The case z ∈ V is similar to z ∈ U . Therefore, (55) is proved. Next, we show that the third term in the RHS of (54) also vanishes, since This equality follows trivially by the fact that for any z U * and . Next, we show that the fourth term in the RHS of (54) satisfies Let . First, let us show that in the unmodified tree, . If x = 0, then this holds by definition of phs. Next, suppose that x ≠ 0. Since , Proposition 1 implies that . Hence also in this case. Next, to analyze in the modified tree, let B be the set of ancestors of u which are in 𝒯 u . By (47) of Lemma 10 , these ancestors b ∈ B are the only inner nodes whose states at location I get modified. Specifically, they satisfy . Hence, if , and otherwise. Let us calculate the number of pairs whose LCA is in B . Denote , and let , where b l +1 is the parent of b l , and b 1 is the parent of u . Observe that b l has 2 l −1 leaves in each of its branches. However, since u is in one of the branches and , there are only 2 l −1 relevant leaves in one of the branches. The total number of pairs is thus This proves (57). Similarly, for the fifth term in the RHS of (54), Next, we show that the last term in the RHS of (54) vanishes, This follows from s u = x ≠ y = s v , which implies that phs( u, v ) = phs ′ ( u, v ) = 0. Finally, we show that the sixth term in the RHS of (54) satisfies Observe that Let z ∈ U * . By Proposition 1 , s LCA( z,u ) = x ≠0, and thus phs( z, u ) = 0. In the modified tree, , and thus also phs ′ ( z, u ) = 0. Next, since s v ≠ s z , we have also phs( z, v ) = 0. In the modified tree, . In addition, since z and v are from different sides of the tree, their LCA is the root which is unmutated. Hence, phs ′ ( z, v ) = 1. Now, let z ∈ V * . By a similar argument, it follows that phs( z, v ) = phs ′ ( z, v ) = phs( z, u ) = 0 and phs ′ ( z, u ) = 1. Hence, (61) reads from which (60) follows. Plugging (55), (56), (57), (58) and (59) and (60) into (54) proves the lemma. I.6 Proof of Lemma 9 (probability distribution of E i ( u )) Proof of Lemma 9 . Let . If , then e u = 0 by definition. Next, consider the case . By Bayes’ law, By definition, A i ( u ) is e u edges above u . Let b be the parent of A i ( u ). Then b is unmutated and A i ( u ) is mutated. Since the depth of b is d − e u − 1, we have In addition, The lemma follows by combining the above three equations. J Generation of Random Ground-Truth Tree Topologies We generated ground-truth trees using the birth-death process implemented in the Cassiopeia package [ JKQ + 20 ], specifying parameters for birth and death rates to model realistic lineage structures observed in single-cell tracing experiments. In our simulations, the time between cell divisions follows an exponential distribution with a birth rate initialized to 2. Tree lineages die at times drawn from an exponential distribution with a rate fixed at 0.75. Upon a cell division, its fitness changes with probability 50%. If the fitness changes, then the birth rate is multiplied by 1.1 z where z is drawn from a Normal distribution with mean 0.5 and standard deviation 0.25. The tree grows until the cell population reaches a predefined value, and then it terminates. The times are rescaled such that the end of the experiment is set to have time τ = 1. Finally, a set of n leaves is subsampled. This subset induces an underlying ground-truth topology, denoted by 𝒯 GT . K Additional Simulation Results As described in Section 4 , our PHS-based test statistic requires the knowledge of the model parameter q . In our simulations, we assumed that this parameter is known for reasons explained in Appendix A.1. Moreover, empirically, our test statistic is not sensitive to the exact value of q . To illustrate this point, we calculated cPHS with three alternatives for the value of q : (i) the actual q ; (ii) a noisy version of it, q + 𝒩 (0, 0.3 q ), trimmed to the range [0.1 q , 1]; and (iii) q = 1 /m , which corresponds to uniform probability distribution of mutation states, namely q i = 1 /m for i = 1, …, m . The result, under the same setting of Figure 3 , is depicted in Figure 10 . As shown, the performance of cPHS is only slightly affected by an inaccurate estimation of q . Download figure Open in new tab Figure 10: AUC values for cPHS with the true value of q , a noisy value of it, and with q = 1 /m that corresponds to uniform mutation probabilities ( q i = 1 /m ). The setting is the same as in Figure 3 (right). In the main text, Figures 2 , 3 , and 4 summarize the outcomes obtained from 7000 reconstructed trees evaluated under each parameter configuration. Figure 11 presents histograms that distinguish the counts of accurately versus inaccurately reconstructed trees within this sample of 7000. In the upper two panels of Figure 11 , the distribution of accurate and inaccurate reconstructions is approximately balanced. In contrast, the lower panel exhibits a pronounced imbalance: the number of inaccurate reconstructions markedly exceeds that of accurate ones. This is due to the fact that the cutoff value ϵ in the lower panel is smaller. Download figure Open in new tab Figure 11: Histograms of accurate and inaccurate reconstructed trees, according to the hypothesis testing (3), as a function of sequence length k and number of mutated states m . Upper left panel: d = d RF and ϵ = 2 / 3. Upper right panel: d = d P and ϵ = 0.01. Lower panel: d = d RF and ϵ = 1 / 2. Next, Figures 12 , 13 , 14 , 15 , and 16 complement Fig. 1 of the main text by showing results for additional cutoff values ϵ and distance functions. Specifically, Figures 12 and 13 show the result for the same parsimony distance function ( d = d P ), but with the cutoff values ϵ = 0.001 and ϵ = 0.03, respectively. Figures 14 , 15 , and 16 show results for the normalized RF ( d = d RF ), likelihood ( d = d L ) and triplets ( d = d tri ) distance functions, with the cutoff values ϵ = 2 / 3, ϵ = 0.01, and ϵ = 1 / 2, respectively. In all these results, the conclusion of Fig. 1 holds: only cPHS exhibits a good separation between the distribution of accurate and inaccurate trees. Download figure Open in new tab Figure 12: Kernel density estimates for the four accuracy measures, with the same setting as in Fig. 1 , but with ϵ = 0.001. Download figure Open in new tab Figure 13: Kernel density estimates for the four accuracy measures, with the same setting as in Fig. 1 , but with ϵ = 0.03. Download figure Open in new tab Figure 14: Kernel density estimates for the four accuracy measures, with the same setting as in Fig. 1 , but with the normalized RF as a distance function ( d = d RF ) and ϵ = 2 / 3. Download figure Open in new tab Figure 15: Kernel density estimates for the four accuracy measures, with the same setting as in Fig. 1 , but with the likelihood as a distance function ( d = d L ) and ϵ = 0.01. Download figure Open in new tab Figure 16: Kernel density estimates for the four accuracy measures, with the same setting as in Fig. 1 , but with the triplets as a distance function ( d = d tri ) and ϵ = 1 / 2. Figures 17 and 18 complement Figures 2 and 3 of the main text, by showing results for the triplets and likelihood distance functions, with cutoff values ϵ = 1 / 5 and ϵ = 0.01, respectively. Figures 19 , 20 , 21 and 22 complement Figures 3 and 4 of the main text but showing results for an additional cutoff values ϵ and additional ranges of parameters. Specifically, Figures 19 and 20 show results for the same normalized RF function ( d = d RF ), but with the cutoff value ϵ = 2 / 3. In this case, the fixed threshold of cPHS is set to t = 10 −11 (as opposed to t = 10 −3 in the main text) to accommodate the higher cutoff value ϵ , thereby permitting a greater number of reconstructed trees to be accepted. Finally, Figures 21 and 22 show results across ranges of two other model parameters: the mutation probability at a leaf ρ , and the number of observed leaves n . Across all these settings, the qualitative conclusions remain consistent with those presented in the main text. Download figure Open in new tab Figure 17: ROC (top) and Precision-Recall (bottom) curves with the same setting as in Fig. 2 except for the distance function. Left: triplets distance ( d = d tri ) with ϵ = 1 / 5 (left panel). Right: likelihood distance ( d = d L ) with ϵ = 0.01. Download figure Open in new tab Figure 18: AUC values with the same setting as in Fig. 3 , except for the distance function. Left: triplets distance ( d = d tri ) with ϵ = 1 / 5. Right: likelihood distance ( d = d L ) with ϵ = 0.01. Download figure Open in new tab Figure 19: AUC values with the same setting as in Fig. 3 , except for the cutoff value which is ϵ = 2 / 3. Download figure Open in new tab Figure 20: Balanced accuracy of the test statistics with the same setting as in Fig. 4 , except for the cutoff value which is ϵ = 2 / 3. Download figure Open in new tab Figure 21: AUC values with the same setting as in Fig. 3 , as a function of ρ (left) and n (right). Download figure Open in new tab Figure 22: Balanced accuracy of the test statistics with the same setting as in Fig. 4 , as a function of ρ (left) and n (right). Funder Information Declared Israel Science Foundation , 2362/22 European Union, https://ror.org/019w4f821 , 101089213 Footnotes { pini.zilber{at}weizmann.ac.il , boaz.nadler{at}weizmann.ac.il } sprillo{at}berkeley.edu nir.yosef{at}weizmann.ac.il ↵ 1 A Python implementation of our method will be made publicly available upon publication. References [Alb05]. ↵ Victor A Albert . Parsimony, phylogeny, and genomics . OUP Oxford , 2005 . [BH95]. ↵ Yoav Benjamini and Yosef Hochberg . Controlling the false discovery rate: a practical and powerful approach to multiple testing . Journal of the Royal statistical society: series B (Methodological) , 57 ( 1 ): 289 – 300 , 1995 . OpenUrl CrossRef PubMed Web of Science [BHY09]. ↵ Yoav Benjamini , Ruth Heller , and Daniel Yekutieli . Selective inference in complex research . Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences , 367 ( 1906 ): 4255 – 4271 , 2009 . OpenUrl CrossRef PubMed [CS65]. ↵ Joseph H Camin and Robert R Sokal . A method for deducing branching sequences in phylogeny . Evolution , pages 311 – 326 , 1965 . [DJS86]. ↵ William HE Day , David S Johnson , and David Sankoff . The computational complexity of inferring rooted phylogenies by parsimony . Mathematical biosciences , 81 ( 1 ): 33 – 42 , 1986 . OpenUrl CrossRef Web of Science [ESSW99]. ↵ Péter L Erdős , Michael A Steel , László A Székely , and Tandy J Warnow . A few logs suffice to build (almost) all trees (i) . Random Structures & Algorithms , 14 ( 2 ): 153 – 184 , 1999 . OpenUrl [FDM+21]. ↵ Jean Feng , William S DeWitt , Aaron McKenna , Noah Simon , Amy D Willis , and Frederick A Matsen IV . Estimation of cell lineage trees by maximum-likelihood phylogenetics . The annals of applied statistics , 15 ( 1 ): 343 , 2021 . OpenUrl PubMed [Fel81]. ↵ Joseph Felsenstein . Evolutionary trees from dna sequences: a maximum likelihood approach . Journal of molecular evolution , 17 : 368 – 376 , 1981 . OpenUrl CrossRef PubMed Web of Science [Fel04]. ↵ Joseph Felsenstein . Inferring Phylogenies . Sinauer Associates . Sunderland, Massachusetts , 2004 . [Fit71]. ↵ Walter M Fitch . Toward defining the course of evolution: minimum change for a specific tree topology . Systematic Biology , 20 ( 4 ): 406 – 416 , 1971 . OpenUrl CrossRef PubMed [Fit77]. ↵ Walter M Fitch . On the problem of discovering the most parsimonious tree . The American Naturalist , 111 ( 978 ): 223 – 257 , 1977 . OpenUrl CrossRef Web of Science [GGH+21]. ↵ Wuming Gong , Alejandro A Granados , Jingyuan Hu , Matthew G Jones , Ofir Raz , Irepan Salvador-Martínez , Hanrui Zhang , Ke-Huan K Chow , Il-Youp Kwak , Renata Retkute , et al. Benchmarked approaches for reconstruction of in vitro cell lineages and in silico models of c. elegans and m. musculus developmental trees . Cell systems , 12 ( 8 ): 810 – 826 , 2021 . OpenUrl PubMed [Gus97]. ↵ Dan Gusfield . Algorithms on Strings, Trees, and Sequences: Computer Science and Computational Biology . Cambridge University Press , 1997 . [Hal04]. ↵ Barry G Hall . Phylogenetic trees made easy . Sinauer Associates . Sunderland, Massachusetts , 2004 . [Hen99]. ↵ Willi Hennig . Phylogenetic systematics . University of Illinois Press , 1999 . [JKQ+20]. ↵ Matthew G Jones , Alex Khodaverdian , Jeffrey J Quinn , Michelle M Chan , Jeffrey A Hussmann , Robert Wang , Chenling Xu , Jonathan S Weissman , and Nir Yosef . Inference of single-cell phylogenies from lineage tracing data using cassiopeia . Genome biology , 21 : 1 – 27 , 2020 . OpenUrl CrossRef PubMed [KVO18]. ↵ Lennart Kester and Alexander Van Oudenaarden . Single-cell transcriptomics meets lineage tracing . Cell stem cell , 23 ( 2 ): 166 – 179 , 2018 . OpenUrl CrossRef PubMed [MG19]. ↵ Aaron McKenna and James A Gagnon . Recording development with single cell dynamic lineage tracing . Development , 146 ( 12 ): dev169730 , 2019 . OpenUrl Abstract / FREE Full Text [NK00]. ↵ Masatoshi Nei and Sudhir Kumar . Molecular evolution and phylogenetics . Oxford University Press , 2000 . [PAW+25]. ↵ Sebastian Prillo , Kevin An , Wilson Wu , Ivan Kristanto , Matthew G Jones , Yun S Song , and Nir Yosef . Tree reconstruction guarantees from crispr-cas9 lineage tracing data using neighbor-joining . In International Conference on Research in Computational Molecular Biology , pages 376 – 380 . Springer , 2025 . [PRYS23]. ↵ Sebastian Prillo , Akshay Ravoor , Nir Yosef , and Yun S Song . ConvexML: Scalable and accurate inference of single-cell chronograms from crispr/cas9 lineage tracing data . biorxiv , 2023 . [RF81]. ↵ David F Robinson and Leslie R Foulds . Comparison of phylogenetic trees . Mathematical biosciences , 53 ( 1-2 ): 131 – 147 , 1981 . OpenUrl CrossRef PubMed Web of Science [San75]. ↵ David Sankoff . Minimal mutation trees of sequences . SIAM Journal on Applied Mathematics , 28 ( 1 ): 35 – 42 , 1975 . OpenUrl CrossRef PubMed Web of Science [SL19]. ↵ Ken Sugino and Tzumin Lee . Robust reconstruction of crispr and tumor lineage using depth metrics . bioRxiv : 609107 , 2019 . [SN87]. ↵ Naruya Saitou and Masatoshi Nei . The neighbor-joining method: a new method for reconstructing phylogenetic trees . Molecular biology and evolution , 4 ( 4 ): 406 – 425 , 1987 . OpenUrl CrossRef PubMed Web of Science [SR06]. ↵ Sagi Snir and Satish Rao . Using max cut to enhance rooted trees consistency . IEEE/ACM transactions on computational biology and bioinformatics , 3 ( 4 ): 323 – 333 , 2006 . OpenUrl [SS22]. ↵ Sophie Seidel and Tanja Stadler . Tidetree: a bayesian phylogenetic framework to estimate single-cell trees and population dynamic parameters from genetic lineage tracing data . Proceedings of the Royal Society B , 289 ( 1986 ): 20221844 , 2022 . OpenUrl CrossRef PubMed [SSCR23]. ↵ Palash Sashittal , Henri Schmidt , Michelle Chan , and Benjamin J Raphael . Startle: a star homoplasy approach for crispr-cas9 lineage tracing . Cell Systems , 14 ( 12 ): 1113 – 1121 , 2023 . OpenUrl PubMed [War18]. ↵ Tandy Warnow . Computational phylogenetics: an introduction to designing methods for phylogeny estimation . Cambridge University Press , 2018 . [Wil19]. ↵ Daniel J Wilson . The harmonic mean p-value for combining dependent tests . Proceedings of the National Academy of Sciences , 116 ( 4 ): 1195 – 1200 , 2019 . OpenUrl Abstract / FREE Full Text [WL11]. ↵ Edward Orlando Wiley and Bruce S Lieberman . Phylogenetics: theory and practice of phylogenetic systematics . John Wiley & Sons . Hoboken, New Jersey , 2011 . [WZKY23]. ↵ Robert Wang , Richard Zhang , Alex Khodaverdian , and Nir Yosef . Theoretical guarantees for phylogeny inference from single-cell lineage tracing . Proceedings of the National Academy of Sciences , 120 ( 12 ): e2203352120 , 2023 . OpenUrl CrossRef PubMed [YJN+22]. ↵ Dian Yang , Matthew G Jones , Santiago Naranjo , William M Rideout , Kyung Hoi Joseph Min , Raymond Ho , Wei Wu , Joseph M Replogle , Jennifer L Page , Jeffrey J Quinn , et al. Lineage tracing reveals the phylodynamics, plasticity, and paths of tumor evolution . Cell , 185 ( 11 ): 1905 – 1923 , 2022 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted July 30, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Can you trust your reconstructed lineage tree? A homoplasy-based approach for irreversible evolution Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Can you trust your reconstructed lineage tree? A homoplasy-based approach for irreversible evolution Pini Zilber , Sebastian Prillo , Nir Yosef , Boaz Nadler bioRxiv 2025.07.27.667007; doi: https://doi.org/10.1101/2025.07.27.667007 Share This Article: Copy Citation Tools Can you trust your reconstructed lineage tree? A homoplasy-based approach for irreversible evolution Pini Zilber , Sebastian Prillo , Nir Yosef , Boaz Nadler bioRxiv 2025.07.27.667007; doi: https://doi.org/10.1101/2025.07.27.667007 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7642) Biochemistry (17715) Bioengineering (13907) Bioinformatics (42005) Biophysics (21472) Cancer Biology (18624) Cell Biology (25534) Clinical Trials (138) Developmental Biology (13390) Ecology (19935) Epidemiology (2067) Evolutionary Biology (24356) Genetics (15617) Genomics (22529) Immunology (17753) Microbiology (40437) Molecular Biology (17200) Neuroscience (88697) Paleontology (667) Pathology (2840) Pharmacology and Toxicology (4829) Physiology (7653) Plant Biology (15171) Scientific Communication and Education (2046) Synthetic Biology (4304) Systems Biology (9827) Zoology (2272)