Critical threshold for ancestral reconstruction by maximum parsimony on general phylogenies ...

24 Pages·2017·0.25 MB·English

Critical threshold for ancestral reconstruction
by maximum parsimony on general
7
1
0 phylogenies ∗
2

g
u Sebastien Roch† Kun-Chieh Wang‡
A

8 August 9, 2017

]
R
P
.
h
Abstract
t
a
m
We consider the problem of inferring an ancestral state from observa-
[ tions at the leaves of a tree, assuming the state evolves along the tree ac-

1 cording to a two-state symmetric Markov process. We establish a general
v branching rate condition under which maximum parsimony, a common re-
4
constructionmethodrequiringonlytheknowledgeofthetree,succeedsbet-
2
5 ter than random guessing uniformly in the depth of the tree. We thereby
2 generalize previous results of [GS10, ZSYL10]. Our results apply to both
0
deterministicandi.i.d.edgeweights.
.
8
0
7
1
:
v
i
X
r
a
∗Keywords: Markov model on a tree, reconstruction problem, branching number, parsimony
principle,applicationtophylogenetics.
†DepartmentofMathematics,UW–Madison. WorksupportedbyNSFgrantsDMS-1248176,
DMS-1149312(CAREER)andDMS-1614242.
‡Department of Mathematics, UW–Madison. Supported by NSF grant DMS-1248176 to
S.Roch.
1 Introduction
Ancestral reconstruction In biology, the inferred evolutionary history of or-
ganismsisdepictedbyaphylogenetictree,thatis,arootedtreewhosebranchings
indicatepastspeciationeventswiththeleavesrepresentinglivingspecies[Fel04].
Theevolutionoffeatures,suchasthenucleotideatagivenpositioninthegenome
ofaspecies,thepresenceorabsenceofaproteinorthenumberofhornsinalizard,
iscommonlyassumedtofollowMarkoviandynamicsalongthetree[Ste16]. That
is, on each edge, the state of the feature changes according to a Markov process;
at bifurcations, two independent copies of the feature evolve along the outgoing
edgesstartingfromthestateatthebranchingpoint.
Hereweconsidertheproblemofinferringanancestralstatefromobservations
of a feature at the leaves of a known phylogenetic tree. We refer to this problem,
which has important applications in biology [Tho04, Lib07, Mos04, DMR11], as
theancestralreconstructionproblem. Manyrigorousresultshavebeenobtainedin
previous work; see, e.g., [KS66, BRZ95, Iof96, TS97, Mos98, EKPS00, Mos01,
MP03, BCMR06, Sly09, GS10, ZSYL10, BST10, BVVW11, Sly11, FRa, FRb]
forapartiallist. Typically,oneseeksanestimatoroftherootstatewhichisstrictly
superiortorandomguessing—uniformlyinthedepthofthetree—underauniform
prior on the root [Mos01]. Whether such an estimator exists has been shown to
depend on a trade-off between the mixing rate of the Markov process and the
growthrateofthetree. Insomecases,forinstanceintwo-statesymmetricMarkov
chains on binary trees [KS66, Iof96] as well as on more general trees [EKPS00],
sharpthresholdshavebeenestablished. Wefocusonthatspecialcasehere.
The threshold of maximum parsimony Maximum parsimony is an ancestral
reconstructionmethodwithalonghistoryinevolutionarybiology. Seee.g.[Fel04,
Ste16]. Itsunderlyingprincipleissimpleandintuitive: itcallsforassigningastate
to each internal vertex in such a way as to minimize the total number of changes
along the edges; the resulting state at the root (possibly not unique) is the desired
reconstructedstate. Anadvantageofmaximumparsimony,whichaccountspartly
for its use in practice, is that it only requires knowledge of the tree—not of the
substitution probabilities along the edges. The latter can be difficult to estimate
accurately from data. In fact maximum parsimony is equivalent to maximum
likelihood under those circumstances [TS97]. Another root state estimator with
thatpropertyisrecursivemajority,studiedin[Mos98,Mos04].
From its very definition, one might expect maximum parsimony to perform
wellwhentheprobabilitiesofsubstitutionalongtheedgesare“sufficientlysmall.”
1
This statement was confirmed rigorously for the two-state symmetric model on
complete binary trees with constant mutation probability in [ZSYL10] (as first
conjectured in [Ste89]). Below a critical probability, maximum parsimony does
indeed beat random guessing. A related result was obtained in [GS10] under
a common random tree model known as the Yule tree. Here we substantially
generalize both of these results. We give a general condition on the branching
rate under which maximum parsimony succeeds. Our condition is related to the
branching number, which plays an important in many processes on trees. See
e.g.[LP16].
Definitions In order to state our results formally, we begin with some defini-
tions. Let T = (V,E) be an infinite complete binary tree rooted at ρ. That is, all
vertices of T have exactly two children; in particular, T has no leaf. Every edge
e ∈ E isassignedaweightθ ∈ [0,1]. Ife = (x,y)wherexistheparentofy,we
e
alsowriteθ = θ . Weusethenotationx ≤ y toindicatethatxisanancestorofy
y e
andwewritex < y forx ≤ y andx (cid:54)= y. Wealsolets(x)bethesiblingofx (cid:54)= ρ.
UndertheCavender-Farris(CF)model[Cav78,Far73]overT andθ = (θ ) ,
y y∈V
alsoknownasNeyman2-statemodel[Ney71],weassociatetoeachvertexx ∈ V
a state σ ∈ {0,1} as follows. The state at the root, σ , is picked uniformly
x ρ
at random in {0,1}. Recursively, if y has parent x, state σ is equal to σ with
y x
probabilityθ ,otherwiseitispickeduniformlyatrandomin{0,1}. Welet
y
1−θ
p = P [σ (cid:54)= σ ] = y,
y T,θ x y
2
be the probability of a substitution on edge e = (x,y). Here P stands for the
T,θ
probabilityoperatorundertheCFmodeloverT andθ. TheCFmodelisequivalent
to the Fortuin-Kastelyn random cluster representation of the ferromagnetic Ising
modelonT withafreeboundary. Seee.g.[EKPS00]andreferencestherein.
Informally, a root state estimator is a map returning a (possibly randomized)
guessforthestateoftheroot,basedontheknowledgeofthestatesonanobserved
cutset.
Definition 1 (Cutset). A cutset is a minimal, finite set of vertices π ⊆ V such
that all infinite self-avoiding paths starting at ρ must visit π. We let C(T ) be the
collectionofallcutsetsofT . WedenotebyT π = (Vπ,Eπ)thefinitetreeobtained
fromT afterremovingalldescendantsoftheverticesinπ.
As mentioned above, our focus in this work is on a root state estimator known in
phylogenetics as maximum parsimony. The parsimony principle dictates that one
2
assignstoeachvertexx(ancestortotheobservedcutsetπ)astateσˆ suchthatthe
x
overallnumberofchangesalongtheedgesofT π,namely,
(cid:88)
1{σˆ (cid:54)= σˆ },
x y
(x,y)∈Eπ
is minimized, where we let by default σˆ = σ for all z ∈ π. In case both 0 and
z z
1 can be obtained in this way as root state, a uniformly random value in {0,1} is
returned. We let RAπ be the reconstruction accuracy of maximum parsimony,
T,θ
i.e.theprobabilitythatitcorrectlyreconstructstherootstate.
Main result: deterministic weights In our main result, we give conditions
under which the reconstruction accuracy of maximum parsimony is uniformly
boundedawayfrom1/2.
Theorem 1 (Reconstruction accuracy of parsimony: deterministic weights). Let
T = (V,E) be an infinite complete binary tree with edge weights θ = (θ )
z z∈V
satisfying
inf θ > 0, (1)
z
z∈V
and
(cid:40) (cid:41)
(cid:88) (cid:89) 3
sup κ > 0 : inf κ−1θ > 0 > . (2)
z
π∈C(T) 2
x∈π ρ(cid:54)=z≤x
Then the reconstruction accuracy of maximum parsimony on T is bounded away
from1/2,thatis,
1
inf RAπ > .
π∈C(T) T,θ 2
Although T is an infinite complete binary tree, since we allow arbitrary cutsets,
our results in fact apply implicitly to general finite trees. Condition (2) is closely
related to the branching number, a generalized notion of branching rate which
plays a key role in the analysis of many stochastic processes on trees and tree-
like graphs. See e.g. [Per99, LP16] for a survey; see also [Lyo89, Lyo90]. The
followingexampleprovidessomeintuitioninaspecialcase.
Example1(Fixededgeweights). Asasimple illustration,observethat,whenall
weights are equal to θ ∈ (0,1], the supremum in (2) is attained for κ = 2θ .
∗ ∗
Indeed,thesumin(2)whenκ = 2θ simplifiesto
∗
(cid:88) (cid:89) (cid:88) (cid:89) 1
(2θ )−1θ = = 1 > 0
∗ z
2
x∈π ρ(cid:54)=z≤x x∈π ρ(cid:54)=z≤x
3
for any cutset π, where the equality can be proved by induction on the graph
distancefromtheroottothefurthestvertexinπ. Ontheotherhand,lettingπ be
n
thecutsetofallverticesatgraphdistancenfromρ,foranyε > 0itholdsthat
(cid:88) (cid:89)
((2+ε)θ )−1θ = (2+ε)−n ·2n → 0,
∗ z
x∈πnρ(cid:54)=z≤x
as n → +∞. Hence, in this case, condition (2) reduces to 2θ > 3/2, that is,
∗
θ > 3/4. Intermsofsubstitutionprobability,thisisp = 1−θ∗ < 1/8.
∗ ∗ 2
Theargumentintheexampleaboveleadstothefollowingcorollary.
Corollary 1 (All substitution probabilities below the threshold). Let the substitu-
tion probabilities (p ) satisfy p ∈ [0,p ] for all z ∈ V, for some p < 1/8.
z z∈V z ∗ ∗
Then the reconstruction accuracy of maximum parsimony on T is bounded away
from1/2.
In general, Theorem 1 cannot be improved in the following sense. It was shown
in [ZSYL10, Theorem 4.1] that when p = p > 1/8, that is, θ = 1 − 2p =
z ∗ z z
θ < 3/4forallz then
∗
1
infRAπn = ,
n T,θ 2
where π is defined in Example 1. On the other hand, it is not known whether
n
thereconstructionaccuracynecessarilyconvergesto1/2if(2)isnotsatisfied. We
leavethisasanopenproblem.
Zhang et al. [ZSYL10] also established the special case of Theorem 1 when
p = p < 1/8 for all z and π = π . Their proof proceeds through a careful
z ∗ n
analysis of the limit of a recurrence for RAπn first derived in [Ste89]. Our more
T,θ
general result follows from a softer argument which relies on the instability of a
fixedpointofthisrecurrencecorrespondingtoasymptoticreconstructionaccuracy
1/2. A more detailed proof sketch is given in Section 2 following some prelim-
inaries. We note that our proof method may be of more general interest, e.g., to
extend the results beyond the two-state case where the higher dimensionality of
the system may complicate significantly the derivation of an explicit limit even
whenedgeweightsareconstant.
Main result: i.i.d. weights We also obtain a related result in the case of edge
weights that are i.i.d. No lower bound on the weights is needed in this case (and,
hence,itdoesnotimmediatelyfollowfromTheorem1).
4
Theorem 2 (Reconstruction accuracy of parsimony: i.i.d. weights). Let T =
(V,E) be an infinite complete binary tree with edge weights θ = (θ ) drawn
z z∈V
i.i.d. from a distribution Θ over (0,1]. Let µ be the mean of Θ and assume that
µ > 3/4. Then,foranyδ > 0,thereisε > 0suchthat
1
inf RAπ > (1+ε),
π∈C(T) T,θ 2
withprobabilityatleast1−δ.
The previous theorem covers in particular the case of the pure birth process, or
Yuletree,whichisapopularrandommodelofphylogenetictrees. Seee.g.[SS03,
Ste16]. Inthatcase,θ = e−2Tz,whereT isanexponentialwithrateλ. Toderive
z z
thecorrespondingthreshold,wenotethat
(cid:90) ∞ λ
µ = E[θ ] = e−2tλe−λtdt = .
z
λ+2
0
Then
3
µ > ⇐⇒ λ > 6,
4
whichisconsistentwiththeresultsof[GS10].
2 Preliminaries
Our proofs are based on a recurrence for the reconstruction accuracy. After re-
viewingthisrecurrence,wegiveaproofsketchofTheorem1.
2.1 Fitch method
The maximum parsimony solution can be computed efficiently by dynamic pro-
gramming,whichinthiscontextisreferredtoastheFitchmethod.
Definition 2 (Fitch method [Fit71, Har73]). The Fitch method recursively con-
structs a set S(cid:98) of possible states for each vertex z ∈ Vπ, starting from π, as
z
follows.
(a) Ifz ∈ π,thenS(cid:98)π = {σ }.
z z
5
(b) Ifz ∈/ π andhaschildrenxandy,then
(cid:40)
S(cid:98)π ∩S(cid:98)π, ifS(cid:98)π ∩S(cid:98)π (cid:54)= ∅
S(cid:98)π = x y x y
z S(cid:98)π ∪S(cid:98)π, o.w.
x y
The method returns the maximum parsimony estimator σˆ which is equal to the
ρ
uniquestateinS(cid:98)π if|S(cid:98)π| = 1,andotherwisereturnsauniformlyrandomvaluein
ρ ρ
{0,1}.
What is described above is the bottom-up phase of the Fitch method. A top-
down phase, which we will not require here, then assigns a state to each vertex
in Vπ consistent with a maximum parsimony solution. See e.g. [SS03]. Let π
be an arbitrary cutset on T with states σ , u ∈ π, and let S(cid:98) , z ∈ Vπ, be the
u z
correspondingreconstructedsetsundertheFitchmethod. Wedefine
(cid:104) (cid:105)
απ = P S(cid:98)π = {σ } ,
z T,θ z z
(cid:104) (cid:105)
βπ = P S(cid:98)π = {1−σ } .
z T,θ z z
Under our randomization rule, the reconstruction accuracy of maximum parsi-
monyisgivenby
(cid:104) (cid:105) 1 (cid:104) (cid:105)
RAπ = P S(cid:98)π = {σ } + P S(cid:98)π = {0,1}
T,θ T,θ ρ ρ 2 T,θ ρ
1
= απ + (1−απ −βπ)
ρ 2 ρ ρ
1 1
= + (απ −βπ). (3)
2 2 ρ ρ
2.2 Proof sketch
Fix θ satisfying the assumptions of Theorem 1 and fix a cutset π. To analyze
(απ,βπ),itisnaturaltotakeadvantageoftherecursivenatureofT [Ste89,SS07,
z z
ZSYL10]. Let x and y be the children of z. The event S(cid:98)π = {σ } occurs when
z z
either (a) S(cid:98)π = {σ } and S(cid:98)π = {σ }, or (b) S(cid:98)π = {σ } and S(cid:98)π = {0,1} or
x z y z x z y
vice versa. By the Markov property of the CF model, the random variables S(cid:98)π
x
andS(cid:98)π,whicharefunctionsonlyofthestatesofπ underxandy respectively,are
y
conditionally independent given σ . Hence, letting q = 1−p for u = x,y and
z u u
6
takingintoaccountthepossibilityofamutationalongtheedges(z,x)and(z,y),
itfollowsthat
απ = (q απ +p βπ)(q απ +p βπ)
z x x x x y y y y
+(q απ +p βπ)(1−απ −βπ)
x x x x y y
+(1−απ −βπ)(q απ +p βπ), (4)
x x y y y y
where the first and second lines on the r.h.s. correspond respectively to cases (a)
and(b)above. Similarly,
βπ = (p απ +q βπ)(p απ +q βπ)
z x x x x y y y y
+(p απ +q βπ)(1−απ −βπ)
x x x x y y
+(1−απ −βπ)(p απ +q βπ). (5)
x x y y y y
Thisrecurrencewasfirstderivedin[Ste89,Lemma7.20]. Seealso[Mad95].
Inthecasethatp = pforallu,afixedpointanalysiswasperformedin[Ste89,
u
Theorem 7.22]. It was found that, if p ≥ 1/8, there is a single fixed point
(1/3,1/3) which corresponds informally to “having no information about the
root.” While if p < 1/8, there is an additional fixed point (α∞,β∞) with α∞ >
p p p
β∞. Convergence to (1/3,1/3) in the first case and (α∞,β∞) in the second
p p p
case was established rigorously in [ZSYL10]. One step in [ZSYL10] involved
the derivation of a new recurrence in terms of the quantities απn − βπn and
z z
1−(απn+βπn),whichfacilitatestheanalysisofthelimitinthefixededgeweight
z z
case.
Goingbacktogeneraltrees,asourstartingpointwefurthermodifytherecur-
renceof[ZSYL10]. Forallz ∈ V,wedefine
dπ = απ −βπ, (6)
z z z
uπ = 3(απ +βπ)−2. (7)
z z z
WeshowinProposition1belowthat(dπ,uπ)satisfiesthefollowingrecurrence
z z
(cid:18)4−uπ(cid:19) (cid:18)4−uπ(cid:19)
dπ = y θ dπ + x θ dπ, (8)
z 6 x x 6 y y
3 1
uπ = θ θ dπdπ − uπuπ, (9)
z 2 x y x y 2 x y
for z ∈ Vπ − π with children x,y, as well as the inequalities 0 ≤ dπ ≤ 1 and
z
−1/2 ≤ uπ ≤ 1 for all z ∈ Vπ and the boundary conditions dπ = uπ = 1 for
z z z
7
all z ∈ π. Our choice of parametrization is motivated in part by the fact that the
“no information” fixed point is now at (0,0) and that |dπ|,|uπ| ≤ 1. At a high
z z
level we show that, under the branching rate condition (2), the fixed point (0,0)
is “unstable” and that dπ in particular stays bounded away form 0. That in terms
z
impliesalowerboundonthereconstructionaccuracyas,by(3),wehave
1 1
RAπ = + dπ.
T,θ 2 2 ρ
Thelinkbetweenstabilityandtheweightedbranchingratein(2)iseasilyseen
from(8). First,whenallweightsareequaltoθ andπ = π ,wegetbysymmetry
∗ n
thatcloseto(0,0)tothefirstorder
4
dπ ≈ θ dπ.
z 3 ∗ x
Hence, in that case, the solution can be expected to grow when 4θ > 1, corre-
3 ∗
sponding to the condition derived in Example 1. More generally, we use (8) to
relate the d-value at the root to the d-values on a cutset (see Lemma 1). To deal
with the nonlinear nature of (8) and (9), we control the u-values thanks to the
quadraticformof(9)whichimpliesaquickdecaytowards0(seeLemma3).
2.3 Recurrence
Before proceeding to the proof of our main results, we first establish a basic re-
currence which follows from the work of [ZSYL10]. We give a short proof for
completeness.
Proposition1(Recurrenceandbasicproperties). Thefollowinghold:
- [Boundaryconditions]Forallz ∈ π,
dπ = 1, (10)
z
uπ = 1. (11)
z
- [Recurrence]Forallz ∈ Vπ −π,ifx,y arethechildrenofz,
(cid:18)4−uπ(cid:19) (cid:18)4−uπ(cid:19)
dπ = y θ dπ + x θ dπ, (12)
z 6 x x 6 y y
3 1
uπ = θ θ dπdπ − uπuπ. (13)
z 2 x y x y 2 x y
8
- [Bounds]Forallz ∈ Vπ,
0 ≤ dπ ≤ 1, (14)
z
−1/2 ≤ uπ ≤ 1. (15)
z
Proof. We start with the boundary conditions. By Definition 2(a) and the defini-
tionsofαπ andβπ,wehaveforallz ∈ π
z z
(cid:104) (cid:105)
απ = P S(cid:98)π = {σ } = 1,
z T,θ z z
(cid:104) (cid:105)
βπ = P S(cid:98)π = {1−σ } = 0.
z T,θ z z
Souπ = 3(απ+βπ)−2 = 1anddπ = απ−βπ = 1,whichproves(11)and(10).
z z z z z z
Letz ∈ Vπ−πwithchildrenx,y. WedefineΣ = απ+βπ and∆ = απ−βπ.
u u u u u u
Bythedefinitionsofp andq ,notethat
x x
(cid:18) (cid:19) (cid:18) (cid:19)
1+θ 1−θ 1 1
q απ +p βπ = x απ + x βπ = Σ + θ ∆ ,
x x x x 2 x 2 x 2 x 2 x x
and,similarly,
1 1
p απ +q βπ = Σ − θ ∆ .
x x x x 2 x 2 x x
Hence,by(4)and(5),
1 1 1
απ = (Σ +θ ∆ )(Σ +θ ∆ )+ (Σ +θ ∆ )(1−Σ )+ (1−Σ )(Σ +θ ∆ ),
z 4 x x x y y y 2 x x x y 2 x y y y
and
1 1 1
βπ = (Σ −θ ∆ )(Σ −θ ∆ )+ (Σ −θ ∆ )(1−Σ )+ (1−Σ )(Σ θ ∆ ).
z 4 x x x y y y 2 x x x y 2 x y y y
Addinguptheabovetwoequations,weget
1 1
Σ = Σ Σ + θ θ ∆ ∆ +Σ (1−Σ )+(1−Σ )Σ
z x y x y x y x y x y
2 2
3 1
= Σ +Σ − Σ Σ + θ θ ∆ ∆ . (16)
x y x y x y x y
2 2
Similarly,
1 1
∆ = θ ∆ Σ + θ ∆ Σ +θ ∆ (1−Σ )+θ ∆ (1−Σ )
z x x y y y x x x y y y x
2 2
(cid:18) (cid:19) (cid:18) (cid:19)
1 1
= 1− Σ θ ∆ + 1− Σ θ ∆ ,
y x x x y y
2 2
9